LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
CFAlgoStartTimeCompat.cxx
Go to the documentation of this file.
1 #ifndef RECOTOOL_CFALGOSTARTTIMECOMPAT_CXX
2 #define RECOTOOL_CFALGOSTARTTIMECOMPAT_CXX
3 
5 #include <algorithm>
6 
7 namespace cmtool {
8 
9  //-------------------------------------------------------
11  //-------------------------------------------------------
12  {
14  _w2cm = geou.WireToCm();
15  _t2cm = geou.TimeToCm();
16  SetVerbose(false);
17  SetTimeDist(5.);//in cm
18  }
19 
20  //-----------------------------
22  //-----------------------------
23  {
24  }
25 
26  //----------------------------------------------------------------------------------------------
27  float CFAlgoStartTimeCompat::Float(const std::vector<const cluster::ClusterParamsAlg*> &clusters)
28  //----------------------------------------------------------------------------------------------
29  {
30 
31  // Code-block by Kazu starts
32  // This ensures the algorithm works only if # clusters is > 2 (and not =2)
33  // You may take out this block if you want to allow matching using clusters from only 2 planes.
34  if(clusters.size()==2) return -1;
35  // Code-block by Kazu ends
36 
37  //This algorithm now works for 3 planes: find 3Dstart point from first 2 planes and find
38  //How well that agrees with 3rd plane's start point location.
39 
40  //So first, make sure clusters vector has only 3 elements. If not return -1
41  if ( clusters.size() != 3 )
42  return -1;
43 
45 
46  //Get Start Wire [cm], Time [cm], Plane and Channel for each cluster in the 3 planes.
47  double startWirecm0 = clusters.at(0)->GetParams().start_point.w;
48  double startTimecm0 = clusters.at(0)->GetParams().start_point.t;
49  int startWire0 = int( startWirecm0 / _w2cm );
50  unsigned char Pl0 = clusters.at(0)->GetParams().start_point.plane;
51  UInt_t startChan0 = geo->PlaneWireToChannel(Pl0, startWire0);
52 
53  double startWirecm1 = clusters.at(1)->GetParams().start_point.w;
54  double startTimecm1 = clusters.at(1)->GetParams().start_point.t;
55  int startWire1 = int( startWirecm1 / _w2cm );
56  unsigned char Pl1 = clusters.at(1)->GetParams().start_point.plane;
57  UInt_t startChan1 = geo->PlaneWireToChannel(Pl1, startWire1);
58 
59  double startWirecm2 = clusters.at(2)->GetParams().start_point.w;
60  double startTimecm2 = clusters.at(2)->GetParams().start_point.t;
61  int startWire2 = int( startWirecm2 / _w2cm );
62  unsigned char Pl2 = clusters.at(2)->GetParams().start_point.plane;
63  UInt_t startChan2 = geo->PlaneWireToChannel(Pl2, startWire2);
64 
65  //Get Intersections in pairs:
66  //y and z indicate detector coordinate and numbers indicate planes
67  //used to generate that intersection point
68  double yS01, zS01, yS02, zS02, yS12, zS12;
69 
70  //Check if Channels Intersect. If so get [Y,Z] intersection points
71  bool WireIntersect01 = geo->ChannelsIntersect( startChan0, startChan1, yS01, zS01);
72  bool WireIntersect02 = geo->ChannelsIntersect( startChan0, startChan2, yS02, zS02);
73  bool WireIntersect12 = geo->ChannelsIntersect( startChan1, startChan2, yS12, zS12);
74 
75  // Check for 2-Plane Start-Point compatibility:
76  // If Start times within pre-defined range
77  // AND wires intersect --> Start Point is Good
78  // Measure of compatibility is time-distance
79  double compat01 = 9999, compat02 = 9999, compat12 = 9999;
80  if ( (abs(startTimecm0 - startTimecm1) < _timeDist)
81  and WireIntersect01 )
82  compat01 = 0.01*abs(100*(startTimecm0-startTimecm1));
83  if ( (abs(startTimecm1 - startTimecm2) < _timeDist)
84  and WireIntersect12 )
85  compat12 = 0.01*abs(100*(startTimecm1-startTimecm2));
86  if ( (abs(startTimecm0 - startTimecm2) < _timeDist)
87  and WireIntersect02 )
88  compat02 = 0.01*abs(100*(startTimecm0-startTimecm2));
89 
90  // If no compatibility return -1
91  if ( (compat01 > 1000) and (compat02 > 1000) and (compat12 > 1000) ){
92  if ( _verbose ) { std::cout << "No compatibility among any 2 planes..." << std::endl; }
93  return -1;
94  }
95 
96  // If any of the planes were compatible, return "best"
97  // This is the one where time-separation is smallest
98  int bestcompat;
99  double bestcompatVal=0;
100  if ( (compat01 < compat02) and (compat01 < compat12) and (compat01 < 9999) ){
101  bestcompat = 2;
102  bestcompatVal = compat01;
103  }
104  if ( (compat02 < compat01) and (compat02 < compat12) and (compat02 < 9999) ){
105  bestcompat = 1;
106  bestcompatVal = compat02;
107  }
108  if ( (compat12 < compat02) and (compat12 < compat01) and (compat12 < 9999) ){
109  bestcompat = 0;
110  bestcompatVal = compat12;
111  }
112 
113  //Now we know which two planes have a good match for start point (if any)
114  //Check if the cluster on the 3rd plane is compatible with these two planes
115 
116  // To do this check that reconstructed start point from best
117  // combination is inside cluster 3
118  double minWire3rdClus = 9999, minTime3rdClus = 9999;
119  double maxWire3rdClus = 0, maxTime3rdClus = 0;
120  double hitWireTMP, hitTimeTMP;
121  Double_t Start[3] = {-100., 0., 0.};
122  double StartT3rdPlane = 0;
123  double wireStart3rdPlane = 0;
124  UInt_t Pl3rd = 0;
125  std::vector<util::PxHit> hits;
126 
127  if ( bestcompat == 2 ){
128  StartT3rdPlane = startTimecm2;//(startTimecm0+startTimecm1)/2.;
129  Pl3rd = Pl2;
130  Start[1] = yS01;
131  Start[2] = zS01;
132  hits = clusters.at(2)->GetHitVector();
133  }
134  if ( bestcompat == 0 ){
135  StartT3rdPlane = startTimecm0;//(startTimecm1+startTimecm2)/2.;
136  Pl3rd = Pl0;
137  Start[1] = yS12;
138  Start[2] = zS12;
139  hits = clusters.at(0)->GetHitVector();
140  }
141  if ( bestcompat == 1 ){
142  StartT3rdPlane = startTimecm1;//(startTimecm0+startTimecm2)/2.;
143  Pl3rd = Pl1;
144  Start[1] = yS02;
145  Start[2] = zS02;
146  hits = clusters.at(1)->GetHitVector();
147  }
148 
149  for (auto& h: hits){
150  hitWireTMP = h.w/_w2cm;
151  hitTimeTMP = h.t;
152  if ( hitWireTMP < minWire3rdClus ) { minWire3rdClus = hitWireTMP; }
153  if ( hitWireTMP > maxWire3rdClus ) { maxWire3rdClus = hitWireTMP; }
154  if ( hitTimeTMP < minTime3rdClus ) { minTime3rdClus = hitTimeTMP; }
155  if ( hitTimeTMP > maxTime3rdClus ) { maxTime3rdClus = hitTimeTMP; }
156  }
157  try { wireStart3rdPlane = geo->NearestWire( Start, Pl3rd); }
158  //catch ( ::larutil::LArUtilException &e) {
159  catch ( ::cet::exception &e) {
160  std::cout << e.what() << std::endl;
161  std::cout << "Exception caught!" << std::endl;
162  wireStart3rdPlane = -1;
163  }
164 
165  if (_verbose){
166  std::cout << "Cluster Start Times:" << std::endl;
167  std::cout << "Pl0: " << startTimecm0 << "\tPl1: " << startTimecm1 << "\tPl2: " << startTimecm2 << std::endl;
168  std::cout << "Best Time Separation: " << bestcompatVal << std::endl;
169  std::cout << "Third Plane: " << Pl3rd << std::endl;
170  std::cout << "Time Range for 3rd plane: [ " << minTime3rdClus << ", " << maxTime3rdClus << " ]" << std::endl;
171  std::cout << "Start Time for 3rd Plane: " << StartT3rdPlane << std::endl;
172  std::cout << "Reco z,y start point: [ " << Start[2] << ", " << Start[1] << " ]" << std::endl;
173  std::cout << "Nearest Wire on 3rd Plane: " << wireStart3rdPlane << std::endl;
174  std::cout << "Wire range 3rd Plane: [ " << minWire3rdClus << ", " << maxWire3rdClus << " ]" << std::endl;
175  }
176 
177  //Make sure start earest wire and start time in range of 3rd cluster
178  if ( wireStart3rdPlane == -1 )
179  return -1;
180 
181  if ( (wireStart3rdPlane <= maxWire3rdClus) and (wireStart3rdPlane >= minWire3rdClus) and
182  (StartT3rdPlane <= maxTime3rdClus) and (StartT3rdPlane >= minTime3rdClus) ){
183  if (_verbose) { std::cout << "Match OK. Return: " << bestcompatVal << std::endl; }
184  return 1./bestcompatVal;
185  }
186 
187  if (_verbose) { std::cout << "\n\n***************************\n" << std::endl; }
188 
189  return -1;
190  }
191 
192  //------------------------------
194  //------------------------------
195  {
196 
197  }
198 
199 }
200 #endif
virtual float Float(const std::vector< const cluster::ClusterParamsAlg * > &clusters)
CFAlgoStartTimeCompat()
Default constructor.
Class def header for a class CFAlgoStartTimeCompat.
Double_t TimeToCm() const
bool ChannelsIntersect(raw::ChannelID_t c1, raw::ChannelID_t c2, double &y, double &z) const
Returns an intersection point of two channels.
Double_t WireToCm() const
void hits()
Definition: readHits.C:15
void SetVerbose(bool on)
Function to set verbose output.
geo::WireID::WireID_t NearestWire(geo::Point_t const &point, geo::PlaneID const &planeid) const
Returns the index of wire closest to position in the specified TPC.
virtual void Reset()
Function to reset the algorithm instance, called together with manager&#39;s Reset()
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
Float_t e
Definition: plot.C:34
Namespace collecting geometry-related classes utilities.
void SetTimeDist(double t)
Function to set verbose output.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33