DisambigAlg35t.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // DisambigAlg35t.cxx
4 //
5 // trj@fnal.gov
6 // tjyang@fnal.gov
7 //
8 // description
9 //
10 // Based on Tom Junk's idea of triplets matching
11 //
12 //
13 //
14 ////////////////////////////////////////////////////////////////////////
15 
16 
17 //Framework includes:
18 #include "art/Framework/Principal/Event.h"
19 #include "art/Framework/Principal/Handle.h"
20 #include "art_root_io/TFileService.h"
21 #include "art_root_io/TFileDirectory.h"
22 #include "messagefacility/MessageLogger/MessageLogger.h"
23 
24 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
25 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
26 #include "lardata/Utilities/AssociationUtil.h"
27 #include "lardataobj/RecoBase/Hit.h"
28 #include "lardataobj/RecoBase/Wire.h"
29 #include "larcorealg/Geometry/CryostatGeo.h"
30 #include "larcorealg/Geometry/TPCGeo.h"
31 #include "larcorealg/Geometry/PlaneGeo.h"
32 #include "larcorealg/Geometry/WireGeo.h"
33 #include "DisambigAlg35t.h"
34 #include "larevt/Filters/ChannelFilter.h"
35 
36 #include <map>
37 #include <cmath>
38 #include <vector>
39 #include <iostream>
40 #include <fstream>
41 #include <cstdlib>
42 
43 #include "TH1D.h"
44 
45 namespace dune{
46 
47  DisambigAlg35t::DisambigAlg35t(fhicl::ParameterSet const& pset)
48  : fDBScan(pset.get< fhicl::ParameterSet >("DBScanAlg"))
49  {
50  fTimeCut = pset.get<double>("TimeCut");
51  fDistanceCut = pset.get<double>("DistanceCut");
52  fDistanceCutClu = pset.get<double>("DistanceCutClu");
53  fTimeWiggle = pset.get<double>("TimeWiggle");
54  fColChanWiggle = pset.get<int>("ColChannelWiggle");
55  fDoCleanUpHits = pset.get<bool>("DoCleanUpHits", true);
56  }
57 
58 
59  //----------------------------------------------------------
60  //----------------------------------------------------------
61  void DisambigAlg35t::RunDisambig(detinfo::DetectorClocksData const& clockData,
62  detinfo::DetectorPropertiesData const& detProp,
63  const std::vector< art::Ptr<recob::Hit> > &OrigHits )
64  {
65  fDisambigHits.clear();
66 
67  art::ServiceHandle<geo::Geometry> geo;
68 
69  std::vector<std::vector<art::Ptr<recob::Hit> > > hitsUV(2);
70  std::vector<art::Ptr<recob::Hit> > hitsZ;
71 
72  // TH1D *histu = new TH1D("histu","histu",4000,0,4000);
73  // TH1D *histv = new TH1D("histv","histv",4000,0,4000);
74  // TH1D *histz = new TH1D("histz","histz",4000,0,4000);
75 
76  for (size_t i = 0; i<OrigHits.size(); ++i){
77  switch (OrigHits[i]->View()){
78  case geo::kU:
79  hitsUV[0].push_back(OrigHits[i]);
80  // if (OrigHits[i]->WireID().TPC==1)
81  // histu->Fill(OrigHits[i]->PeakTime()
82  // - detProp.GetXTicksOffset(0,
83  // OrigHits[i]->WireID().TPC,
84  // OrigHits[i]->WireID().Cryostat)
85  // ,OrigHits[i]->Charge());
86  break;
87  case geo::kV:
88  hitsUV[1].push_back(OrigHits[i]);
89  // if (OrigHits[i]->WireID().TPC==1)
90  // histv->Fill(OrigHits[i]->PeakTime()
91  // - detProp.GetXTicksOffset(1,
92  // OrigHits[i]->WireID().TPC,
93  // OrigHits[i]->WireID().Cryostat)
94  // ,OrigHits[i]->Charge());
95  break;
96  case geo::kZ:
97  hitsZ.push_back(OrigHits[i]);
98  // if (OrigHits[i]->WireID().TPC==1)
99  // histz->Fill(OrigHits[i]->PeakTime()
100  // - detProp.GetXTicksOffset(2,
101  // OrigHits[i]->WireID().TPC,
102  // OrigHits[i]->WireID().Cryostat)
103  // ,OrigHits[i]->Charge());
104  break;
105  default:
106  throw cet::exception("DisambigAlg35t")
107  <<": hit view unkonw. \n";
108  }
109  }
110  // std::cout<<histu->GetMean()<<" "<<histv->GetMean()<<" "<<histz->GetMean()<<std::endl;
111  // delete histu;
112  // delete histv;
113  // delete histz;
114  std::vector<std::map<unsigned int, unsigned int> >fHasBeenDisambigedUV(2);
115  const unsigned int ntpc = geo->NTPC();
116  //hits and wireids for DBScan
117  std::vector<std::vector<art::Ptr<recob::Hit> > > allhitsu(ntpc);
118  std::vector<std::vector<art::Ptr<recob::Hit> > > allhitsv(ntpc);
119  std::vector<std::vector<art::Ptr<recob::Hit> > > allhitsz(ntpc);
120  std::vector<std::vector<geo::WireID> > wireidsu(ntpc);
121  std::vector<std::vector<geo::WireID> > wireidsv(ntpc);
122 
123  std::vector<std::vector<unsigned int> > cluidu(ntpc);
124  std::vector<std::vector<unsigned int> > cluidv(ntpc);
125  std::vector<std::vector<unsigned int> > bestwireidu(ntpc);
126  std::vector<std::vector<unsigned int> > bestwireidv(ntpc);
127 
128  //Look for triplets of U,V,Z hits that are common in time
129  //Try all possible wire segments for U and V hits and
130  //see if U, V, Z wires cross
131  for (size_t z = 0; z<hitsZ.size(); ++z){//loop over z hits
132  //find triplets of U,V,Z hits that are common in time
133  unsigned int apaz(0), cryoz(0);
134  fAPAGeo.ChannelToAPA(hitsZ[z]->Channel(), apaz, cryoz);
135 
136  double tz = hitsZ[z]->PeakTime()
137  - detProp.GetXTicksOffset(hitsZ[z]->WireID().Plane,
138  hitsZ[z]->WireID().TPC,
139  hitsZ[z]->WireID().Cryostat);
140  //loop over u hits
141  bool findmatch = false;
142  for (size_t u = 0; u<hitsUV[0].size() && !findmatch; ++u){
143  if (fHasBeenDisambigedUV[0].find(u)!=fHasBeenDisambigedUV[0].end()) continue;
144  unsigned int apau(0), cryou(0);
145  fAPAGeo.ChannelToAPA(hitsUV[0][u]->Channel(), apau, cryou);
146  if (apau!=apaz) continue;
147  double tu = hitsUV[0][u]->PeakTime()
148  - detProp.GetXTicksOffset(0,
149  hitsZ[z]->WireID().TPC,
150  hitsZ[z]->WireID().Cryostat);
151 
152  if (std::abs(tu-tz)<fTimeCut){
153  //find a matched u hit, loop over v hits
154  for (size_t v = 0; v<hitsUV[1].size(); ++v){
155  if (fHasBeenDisambigedUV[1].find(v)!=fHasBeenDisambigedUV[1].end()) continue;
156  unsigned int apav(0), cryov(0);
157  fAPAGeo.ChannelToAPA(hitsUV[1][v]->Channel(), apav, cryov);
158  if (apav!=apaz) continue;
159  double tv = hitsUV[1][v]->PeakTime()
160  - detProp.GetXTicksOffset(1,
161  hitsZ[z]->WireID().TPC,
162  hitsZ[z]->WireID().Cryostat);
163 
164  if (std::abs(tv-tz)<fTimeCut){
165 
166  //find a matched v hit, see if the 3 wire segments cross
167  geo::WireID zwire = geo->ChannelToWire(hitsZ[z]->Channel())[0];
168  std::vector<geo::WireID> uwires = geo->ChannelToWire(hitsUV[0][u]->Channel());
169  std::vector<geo::WireID> vwires = geo->ChannelToWire(hitsUV[1][v]->Channel());
170 
171  unsigned int totalintersections = 0;
172  unsigned int bestu = 0; //index to wires associated with channel
173  unsigned int bestv = 0;
174 
175  for (size_t uw = 0; uw<uwires.size(); ++uw){
176  for (size_t vw = 0; vw<vwires.size(); ++vw){
177  geo::WireIDIntersection widiuv;
178  geo::WireIDIntersection widiuz;
179  geo::WireIDIntersection widivz;
180 
181  if (uwires[uw].TPC!=vwires[vw].TPC) continue;
182  if (uwires[uw].TPC!=zwire.TPC) continue;
183  if (vwires[vw].TPC!=zwire.TPC) continue;
184 
185  if (!geo->WireIDsIntersect(zwire,uwires[uw],widiuz)) continue;
186  if (!geo->WireIDsIntersect(zwire,vwires[vw],widivz)) continue;
187  if (!geo->WireIDsIntersect(uwires[uw],vwires[vw],widiuv)) continue;
188 
189  double dis1 = sqrt(pow(widiuz.y-widivz.y,2)+pow(widiuz.z-widivz.z,2));
190  double dis2 = sqrt(pow(widiuz.y-widiuv.y,2)+pow(widiuz.z-widiuv.z,2));
191  double dis3 = sqrt(pow(widiuv.y-widivz.y,2)+pow(widiuv.z-widivz.z,2));
192  double maxdis = std::max(dis1,dis2);
193  maxdis = std::max(maxdis,dis3);
194 
195  if (maxdis<fDistanceCut){
196  ++totalintersections;
197  bestu = uw;
198  bestv = vw;
199  }
200  }
201  }
202  if (totalintersections==1){
203  allhitsu[uwires[bestu].TPC].push_back(hitsUV[0][u]);
204  allhitsv[vwires[bestv].TPC].push_back(hitsUV[1][v]);
205  allhitsz[hitsZ[z]->WireID().TPC].push_back(hitsZ[z]);
206  wireidsu[uwires[bestu].TPC].push_back(uwires[bestu]);
207  wireidsv[vwires[bestv].TPC].push_back(vwires[bestv]);
208  cluidu[uwires[bestu].TPC].push_back(u);
209  cluidv[vwires[bestv].TPC].push_back(v);
210  bestwireidu[uwires[bestu].TPC].push_back(1+bestu);
211  bestwireidv[vwires[bestv].TPC].push_back(1+bestv);
212  findmatch = true;
213  break;
214  }
215  }//find v hit consistent in time
216  }//loop over all v hits
217  }//find u hit consistent in time
218  }//loop over all u hits
219  }//loop over all z hits
220  //Done finding trivial disambiguated hits
221 
222  if (fDoCleanUpHits){
223  //running DB scan to identify and remove outlier hits
224  // get the ChannelFilter
225  filter::ChannelFilter chanFilt;
226  double CorrectedHitTime = 0;
227  int ChannelNumber = 0;
228  for (size_t i = 0; i<ntpc; ++i){//loop over all TPCs
229  if (!allhitsu[i].size()) continue;
230 
231  //initialize dbscan with all hits in u view
232  fDBScan.InitScan(clockData, detProp, allhitsu[i], chanFilt.SetOfBadChannels(), wireidsu[i]);
233  //run dbscan
234  fDBScan.run_cluster();
235  //fpointId_to_clusterId maps a hit index to a cluster index (which can be negative if the hit is not associated with any clusters)
236  if (allhitsu[i].size()!=fDBScan.fpointId_to_clusterId.size())
237  throw cet::exception("DisambigAlg35t") <<"DBScan hits do not match input hits"<<allhitsu[i].size()<<" "<<fDBScan.fpointId_to_clusterId.size()<<"\n";
238 
239  //Find out which clusters are unique in time
240  std::vector<unsigned int> hitstoaddu;
241  std::vector<unsigned int> dbcluhits(fDBScan.fclusters.size(),0);
242  std::vector< bool > boolVector(fDBScan.fclusters.size(), true);
243  std::map<int, std::pair <double,double> > ClusterStartEndTime;
244  std::map<int, std::pair <int, int> > ClusterStartEndColChan;
245  for(size_t j = 0; j < fDBScan.fpointId_to_clusterId.size(); ++j){
246  // for c2: fDBScan.fpointId_to_clusterId[j]>=0 is always true
247  //if (fDBScan.fpointId_to_clusterId[j]>=0&&fDBScan.fpointId_to_clusterId[j]<fDBScan.fclusters.size()) {
248  if (fDBScan.fpointId_to_clusterId[j]<fDBScan.fclusters.size()) {
249  ++dbcluhits[fDBScan.fpointId_to_clusterId[j]];
250 
251  CorrectedHitTime = allhitsu[i][j]->PeakTime() - detProp.GetXTicksOffset(allhitsu[i][j]->WireID().Plane, allhitsu[i][j]->WireID().TPC, allhitsu[i][j]->WireID().Cryostat);
252  if (ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].first==0 || ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].first > CorrectedHitTime)
253  ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].first = CorrectedHitTime;
254  if (ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].second==0 || ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].second < CorrectedHitTime)
255  ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].second = CorrectedHitTime;
256 
257  ChannelNumber = allhitsz[i][j]->Channel();
258  if (ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].first==0 || ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].first > ChannelNumber)
259  ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].first = ChannelNumber;
260  if (ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].second==0 || ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].second < ChannelNumber)
261  ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].second = ChannelNumber;
262  /*
263  std::cout << "Looking at Cluster " << fDBScan.fpointId_to_clusterId[j] << " in TPC " << (int)i
264  << ", It has start time " << ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].first
265  << ", and end time " << ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].second
266  << ", It has start ColChan " << ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].first
267  << ", and end ColChan " << ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].second
268  << std::endl;
269  //*/
270  }
271  }
272  for ( unsigned int ClusNum = 0; ClusNum < fDBScan.fclusters.size(); ++ClusNum ) {
273  for ( unsigned int ClusIt = 0; ClusIt < fDBScan.fclusters.size(); ++ClusIt ) {
274  if ( dbcluhits[ClusIt] > dbcluhits[ClusNum] ) {
275  // Smaller cluster, so check not in the same time range...
276  if ( ClusterStartEndTime[ClusNum].first > ( ClusterStartEndTime[ClusIt].first + fTimeWiggle )
277  && ClusterStartEndTime[ClusNum].first < ( ClusterStartEndTime[ClusIt].second + fTimeWiggle ) ) {
278  if ( ClusterStartEndTime[ClusNum].second > ( ClusterStartEndTime[ClusIt].first + fTimeWiggle )
279  && ClusterStartEndTime[ClusNum].second < ( ClusterStartEndTime[ClusIt].second + fTimeWiggle ) ) {
280  // Within the same time...now check if same collection channel range...
281  if ( ClusterStartEndColChan[ClusNum].first > ( ClusterStartEndColChan[ClusIt].first + fColChanWiggle )
282  && ClusterStartEndColChan[ClusNum].first < ( ClusterStartEndColChan[ClusIt].second + fColChanWiggle ) ) {
283  if ( ClusterStartEndColChan[ClusNum].second > ( ClusterStartEndColChan[ClusIt].first + fColChanWiggle )
284  && ClusterStartEndColChan[ClusNum].second < ( ClusterStartEndColChan[ClusIt].second + fColChanWiggle ) ) {
285  // Within same time and collection channel range...Bad Cluster?
286  boolVector[ClusNum] = false;
287  }
288  }
289  }
290  }
291  }
292  }
293  /*
294  std::cout << "\nLooking at Cluster Number " << ClusNum << ".\n"
295  << "It had start time " << ClusterStartEndTime[ClusNum].first << ", and end time " << ClusterStartEndTime[ClusNum].second << ".\n"
296  << "It had start ColChan " << ClusterStartEndColChan[ClusNum].first << ", and end ColChan " << ClusterStartEndColChan[ClusNum].second << ".\n"
297  << "The bool vector value for this cluster is..." << boolVector[ClusNum]
298  << std::endl;
299  //*/
300  }
301 
302  //std::cout << "Now going to look at v hits....." << std::endl;
303 
304  //hitstoaddu holds all u hits in the time separated clusters
305  for(size_t j = 0; j < fDBScan.fpointId_to_clusterId.size(); ++j){
306  if (int(fDBScan.fpointId_to_clusterId[j]) == -1 ) continue;
307  if (boolVector[fDBScan.fpointId_to_clusterId[j]] == true ) hitstoaddu.push_back(j);
308  }
309 
310  //Now to do the same for v hits
311  fDBScan.InitScan(clockData, detProp, allhitsv[i], chanFilt.SetOfBadChannels(), wireidsv[i]);
312  fDBScan.run_cluster();
313  if (allhitsv[i].size()!=fDBScan.fpointId_to_clusterId.size())
314  throw cet::exception("DisambigAlg35t") <<"DBScan hits do not match input hits"<<allhitsv[i].size()<<" "<<fDBScan.fpointId_to_clusterId.size()<<"\n";
315 
316  std::vector<unsigned int> hitstoaddv;
317  dbcluhits.clear();
318  boolVector.clear();
319  ClusterStartEndTime.clear();
320  ClusterStartEndColChan.clear();
321  dbcluhits.resize(fDBScan.fclusters.size(),0);
322  boolVector.resize(fDBScan.fclusters.size(),true);
323  for(size_t j = 0; j < fDBScan.fpointId_to_clusterId.size(); ++j){
324  //for c2: fDBScan.fpointId_to_clusterId[j]>=0 is always true
325  //if (fDBScan.fpointId_to_clusterId[j]>=0&&fDBScan.fpointId_to_clusterId[j]<fDBScan.fclusters.size()) {
326  if (fDBScan.fpointId_to_clusterId[j]<fDBScan.fclusters.size()) {
327  ++dbcluhits[fDBScan.fpointId_to_clusterId[j]];
328 
329  CorrectedHitTime = allhitsv[i][j]->PeakTime() - detProp.GetXTicksOffset(allhitsv[i][j]->WireID().Plane, allhitsv[i][j]->WireID().TPC, allhitsv[i][j]->WireID().Cryostat);
330  if (ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].first==0 || ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].first > CorrectedHitTime)
331  ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].first = CorrectedHitTime;
332  if (ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].second==0 || ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].second < CorrectedHitTime)
333  ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].second = CorrectedHitTime;
334 
335  ChannelNumber = allhitsz[i][j]->Channel();
336  if (ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].first==0 || ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].first > ChannelNumber)
337  ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].first = ChannelNumber;
338  if (ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].second==0 || ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].second < ChannelNumber)
339  ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].second = ChannelNumber;
340  /*
341  std::cout << "Looking at Cluster " << fDBScan.fpointId_to_clusterId[j] << " in TPC " << (int)i
342  << ", It has start time " << ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].first
343  << ", and end time " << ClusterStartEndTime[fDBScan.fpointId_to_clusterId[j]].second
344  << ", It has start ColChan " << ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].first
345  << ", and end ColChan " << ClusterStartEndColChan[fDBScan.fpointId_to_clusterId[j]].second
346  << std::endl;
347 
348  //*/
349  }
350  }
351  for ( unsigned int ClusNum = 0; ClusNum < fDBScan.fclusters.size(); ++ClusNum ) {
352  for ( unsigned int ClusIt = 0; ClusIt < fDBScan.fclusters.size(); ++ClusIt ) {
353  if ( dbcluhits[ClusIt] > dbcluhits[ClusNum] ) {
354  // Smaller cluster, so check not in the same time range...
355  if ( ClusterStartEndTime[ClusNum].first > ( ClusterStartEndTime[ClusIt].first + fTimeWiggle )
356  && ClusterStartEndTime[ClusNum].first < ( ClusterStartEndTime[ClusIt].second + fTimeWiggle ) ) {
357  if ( ClusterStartEndTime[ClusNum].second > ( ClusterStartEndTime[ClusIt].first + fTimeWiggle )
358  && ClusterStartEndTime[ClusNum].second < ( ClusterStartEndTime[ClusIt].second + fTimeWiggle ) ) {
359  // Within the same time...now check if same collection channel range...
360  if ( ClusterStartEndColChan[ClusNum].first > ( ClusterStartEndColChan[ClusIt].first + fColChanWiggle )
361  && ClusterStartEndColChan[ClusNum].first < ( ClusterStartEndColChan[ClusIt].second + fColChanWiggle ) ) {
362  if ( ClusterStartEndColChan[ClusNum].second > ( ClusterStartEndColChan[ClusIt].first + fColChanWiggle )
363  && ClusterStartEndColChan[ClusNum].second < ( ClusterStartEndColChan[ClusIt].second + fColChanWiggle ) ) {
364  // Within same time and collection channel range...Bad Cluster?
365  boolVector[ClusNum] = false;
366  }
367  }
368  }
369  }
370  }
371  }
372  /*
373  std::cout << "\nLooking at Cluster Number " << ClusNum << ".\n"
374  << "It had start time " << ClusterStartEndTime[ClusNum].first << ", and end time " << ClusterStartEndTime[ClusNum].second << ".\n"
375  << "It had start ColChan " << ClusterStartEndColChan[ClusNum].first << ", and end ColChan " << ClusterStartEndColChan[ClusNum].second << ".\n"
376  << "The bool vector value for this cluster is..." << boolVector[ClusNum] << std::endl;
377  //*/
378  }
379  //hitstoaddv holds all v hits in the time separated clusters
380  for(size_t j = 0; j < fDBScan.fpointId_to_clusterId.size(); ++j){
381  if (int(fDBScan.fpointId_to_clusterId[j]) == -1 ) continue;
382  if (boolVector[fDBScan.fpointId_to_clusterId[j]] == true ) hitstoaddv.push_back(j);
383  }
384  //std::cout << "Looking back through all hits to see if can add hits? Loops through " << hitstoaddu.size() << " hits for U, and " << hitstoaddv.size() << " for V." << std::endl;
385  //size(allhitsu)=size(allhitsv)
386  for (size_t j = 0; j < allhitsu[i].size(); ++j){
387  bool foundu = false;
388  bool foundv = false;
389  for (size_t k = 0; k<hitstoaddu.size(); ++k){
390  if (j==hitstoaddu[k]) {
391  foundu = true;
392  }
393  }
394  for (size_t k = 0; k<hitstoaddv.size(); ++k){
395  if (j==hitstoaddv[k]) {
396  foundv = true;
397  }
398  }
399  if (foundu&&foundv){
400  //only add each hit once
401  if (fHasBeenDisambigedUV[0].find(cluidu[i][j])==fHasBeenDisambigedUV[0].end()){
402  fDisambigHits.push_back(std::pair<art::Ptr<recob::Hit>, geo::WireID>(allhitsu[i][j],wireidsu[i][j]));
403  fHasBeenDisambigedUV[0][cluidu[i][j]] = bestwireidu[i][j];
404  }
405  if (fHasBeenDisambigedUV[1].find(cluidv[i][j])==fHasBeenDisambigedUV[1].end()){
406 
407  fDisambigHits.push_back(std::pair<art::Ptr<recob::Hit>, geo::WireID>(allhitsv[i][j],wireidsv[i][j]));
408  fHasBeenDisambigedUV[1][cluidv[i][j]] = bestwireidv[i][j];
409  }
410  }
411  }
412  }
413 
414  }//if (fDoCleanUpHits)
415  else{//just take all triplets of hits
416  for (size_t i = 0; i<ntpc; ++i){//loop over all TPCs
417  for (size_t j = 0; j < allhitsu[i].size(); ++j){
418  //only add each hit once
419  if (fHasBeenDisambigedUV[0].find(cluidu[i][j])==fHasBeenDisambigedUV[0].end()){
420  fDisambigHits.push_back(std::pair<art::Ptr<recob::Hit>, geo::WireID>(allhitsu[i][j],wireidsu[i][j]));
421  fHasBeenDisambigedUV[0][cluidu[i][j]] = bestwireidu[i][j];
422  }
423  if (fHasBeenDisambigedUV[1].find(cluidv[i][j])==fHasBeenDisambigedUV[1].end()){
424 
425  fDisambigHits.push_back(std::pair<art::Ptr<recob::Hit>, geo::WireID>(allhitsv[i][j],wireidsv[i][j]));
426  fHasBeenDisambigedUV[1][cluidv[i][j]] = bestwireidv[i][j];
427  }
428  }
429  }
430  }
431 
432  //loop over undisambiguated hits, find the nearest channel of disambiguated hits and determine the correct wire segment.
433  for (size_t i = 0; i<2; ++i){//loop over U and V hits
434  for (size_t hit = 0; hit<hitsUV[i].size(); ++hit){
435  if (fHasBeenDisambigedUV[i].find(hit)!=fHasBeenDisambigedUV[i].end()) continue;
436  unsigned int apa1(0), cryo1(0);
437  fAPAGeo.ChannelToAPA(hitsUV[i][hit]->Channel(), apa1, cryo1);
438  //unsigned int channdiff = 100000;
439  //geo::WireID nearestwire;
440  std::vector<geo::WireID> wires = geo->ChannelToWire(hitsUV[i][hit]->Channel());
441  //std::vector<double> disttoallhits(wires.size(),-1);
442  std::vector<int> nearbyhits(wires.size(),-1);
443  for (auto& u2 : fHasBeenDisambigedUV[i]){
444  unsigned int apa2(0), cryo2(0);
445  fAPAGeo.ChannelToAPA(hitsUV[i][u2.first]->Channel(), apa2, cryo2);
446  if (apa1!=apa2) continue;
447  geo::WireID hitwire = geo->ChannelToWire(hitsUV[i][u2.first]->Channel())[u2.second-1];
448  double wire_pitch = geo->WirePitch(hitwire.Plane,hitwire.TPC,hitwire.Cryostat); //wire pitch in cm
449  for (size_t w = 0; w<wires.size(); ++w){
450  if (wires[w].TPC!= hitwire.TPC) continue;
451  if (nearbyhits[w]<0) nearbyhits[w] = 0;
452  double distance = sqrt(pow((double(hitwire.Wire)-double(wires[w].Wire))*wire_pitch,2)+pow(detProp.ConvertTicksToX(hitsUV[i][u2.first]->PeakTime(),hitwire.Plane,hitwire.TPC,hitwire.Cryostat)-detProp.ConvertTicksToX(hitsUV[i][hit]->PeakTime(),hitwire.Plane,hitwire.TPC,hitwire.Cryostat),2));
453  if (distance<fDistanceCutClu) ++nearbyhits[w];
454  }
455  }
456  unsigned bestwire = 0;
457  int maxnumhits = 0;
458  for (size_t w = 0; w<wires.size(); ++w){
459  if (nearbyhits[w]<0) continue;
460  if (nearbyhits[w]>maxnumhits){
461  bestwire = w;
462  maxnumhits = nearbyhits[w];
463  }
464  }
465  if (maxnumhits){
466  fDisambigHits.push_back(std::pair<art::Ptr<recob::Hit>, geo::WireID>(hitsUV[i][hit],wires[bestwire]));
467  fHasBeenDisambigedUV[i][hit] = 1+bestwire;
468  }
469  else{
470  mf::LogInfo("DisambigAlg35t")<<"Could not find disambiguated hit for "<<std::string(hitsUV[i][hit]->WireID())<<" "<<hitsUV[i][hit]->PeakTime();
471  }
472  }
473  }
474  }
475 
476 } //end namespace apa
ValidateOpDetReco.end
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
Definition: ValidateOpDetReco.py:548
geo::WireIDIntersection::z
double z
z position of intersection
Definition: geo_types.h:805
dune::DisambigAlg35t::fAPAGeo
dune::apa::APAGeometryAlg fAPAGeo
Definition: DisambigAlg35t.h:61
CryostatGeo.h
Encapsulate the construction of a single cyostat.
View
AdcChannelData::View View
Definition: test_AdcChannelData.cxx:23
string
std::string string
Definition: nybbler.cc:12
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
geo::kV
Planes which measure V.
Definition: geo_types.h:130
detinfo::DetectorPropertiesData::GetXTicksOffset
double GetXTicksOffset(int p, int t, int c) const
Definition: DetectorPropertiesData.cc:148
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
vector
struct vector vector
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::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
geo::kZ
Planes which measure Z direction.
Definition: geo_types.h:132
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
dune::DisambigAlg35t::DisambigAlg35t
DisambigAlg35t(fhicl::ParameterSet const &pset)
Definition: DisambigAlg35t.cxx:47
geo::GeometryCore::WirePitch
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:824
dune::DisambigAlg35t::fDisambigHits
std::vector< std::pair< art::Ptr< recob::Hit >, geo::WireID > > fDisambigHits
The final list of hits to pass back to be made.
Definition: DisambigAlg35t.h:55
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
geo::kU
Planes which measure U.
Definition: geo_types.h:129
cluster::DBScanAlg::fpointId_to_clusterId
std::vector< unsigned int > fpointId_to_clusterId
mapping point_id -> clusterId
Definition: DBScanAlg.h:73
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:271
protoana::distance
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
Definition: truepionXsection_module.cc:217
dune::apa::APAGeometryAlg::ChannelToAPA
unsigned int ChannelToAPA(uint32_t chan)
Get number of the APA containing the given channel.
Definition: APAGeometryAlg.cxx:141
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
dune::DisambigAlg35t::RunDisambig
void RunDisambig(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit > > &OrigHits)
Run disambiguation as currently configured.
Definition: DisambigAlg35t.cxx:61
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
WireGeo.h
Encapsulate the geometry of a wire.
detinfo::DetectorPropertiesData::ConvertTicksToX
double ConvertTicksToX(double ticks, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:115
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
geo::GeometryCore::WireIDsIntersect
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, geo::Point_t &intersection) const
Computes the intersection between two wires.
cluster::DBScanAlg::InitScan
void InitScan(const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &allhits, std::set< uint32_t > badChannels, const std::vector< geo::WireID > &wireids=std::vector< geo::WireID >())
Definition: DBScanAlg.cxx:272
Hit.h
Declaration of signal hit object.
PlaneGeo.h
Encapsulate the construction of a single detector plane.
cluster::DBScanAlg::fclusters
std::vector< std::vector< unsigned int > > fclusters
collection of something
Definition: DBScanAlg.h:71
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:283
geo::WireIDIntersection::y
double y
y position of intersection
Definition: geo_types.h:804
Wire.h
Declaration of basic channel signal object.
dune::DisambigAlg35t::fDBScan
cluster::DBScanAlg fDBScan
object that implements the DB scan algorithm
Definition: DisambigAlg35t.h:68
art::get
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
util::quantities::second
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
geo
LArSoft geometry interface.
Definition: ChannelGeo.h:16
filter::ChannelFilter::SetOfBadChannels
std::set< uint32_t > SetOfBadChannels() const
Definition: ChannelFilter.cxx:75
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
TPCGeo.h
Encapsulate the construction of a single detector plane.
Generated by   doxygen 1.8.11