Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DFDCIntersection_factory.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: DFDCIntersection_factory.cc
4 // Created: Tue Oct 30 11:24:53 EDT 2007
5 // Creator: davidl (on Darwin fwing-dhcp95.jlab.org 8.10.1 i386)
6 //
7 
8 #include <cmath>
9 using namespace std;
10 
12 #include "FDC/DFDCGeometry.h"
13 #include "HDGEOMETRY/DGeometry.h"
14 
15 //------------------
16 // init
17 //------------------
19 {
20  MAX_DIST2 = 2.0*2.0;
21  DEBUG_LEVEL=0;
22 
23  // Create skeleton of data structure to hold hits
24  vector<const DFDCHit*> mt_trkhits;
25  vector<vector<const DFDCHit*> > mt_trkhits_by_layer;
26  for(int i=0; i<6; i++)mt_trkhits_by_layer.push_back(mt_trkhits);
27  for(int i=0; i<4; i++)fdchits_by_package.push_back(mt_trkhits_by_layer);
28 
29  return NOERROR;
30 }
31 
32 //------------------
33 // brun
34 //------------------
35 jerror_t DFDCIntersection_factory::brun(JEventLoop *loop, int32_t runnumber)
36 {
37  // Get pointer to DGeometry object
38  DApplication* dapp=dynamic_cast<DApplication*>(eventLoop->GetJApplication());
39  DGeometry *dgeom = dapp->GetDGeometry(runnumber);
40  if (!dgeom->GetFDCWires(fdcwires)){
41  _DBG_<< "FDC geometry not available!" <<endl;
42  USE_FDC=false;
43  }
44 
45  return NOERROR;
46 }
47 
48 //------------------
49 // evnt
50 //------------------
51 jerror_t DFDCIntersection_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
52 {
53  // if (!USE_FDC) return NOERROR;
54 
55  // Clear fdchits_by_package structure
56  for(int package=1; package<=4; package++){
57  for(int layer=1; layer<=6; layer++){
58  fdchits_by_package[package-1][layer-1].clear();
59  }
60  }
61 
62  // Get raw hits
63  vector<const DFDCHit*> fdchits;
64  loop->Get(fdchits);
65 
66  // For events with a very large number of hits, assume
67  // we can't reconstruct them so bail early
68  // Feb. 8, 2008 D.L.
69  if(fdchits.size()>(5+5+1)*24*4){
70  if (DEBUG_LEVEL>0){
71  _DBG_<<"Too many hits in FDC! Intersection point reconstruction in FDC bypassed for event "<<loop->GetJEvent().GetEventNumber()<<endl;
72  }
73  return NOERROR;
74  }
75 
76  // Sort wire hits by package and layer
77  int wire=0,layer=0;
78  for(unsigned int i=0; i<fdchits.size(); i++){
79  const DFDCHit *fdchit = fdchits[i];
80 
81  if(fdchit->type != 0)continue; // filter out cathode hits
82  if(fdchit->gLayer<1 || fdchit->gLayer>24){
83  _DBG_<<"FDC gLayer out of range! ("<<fdchit->gLayer<<" is not between 1 and 24)"<<endl;
84  continue;
85  }
86  // Only include the first hit on a given wire in a given layer
87  if (fdchit->element!=wire || fdchit->gLayer!=layer){
88  int package = (fdchit->gLayer-1)/6 + 1;
89  fdchits_by_package[package-1][(fdchit->gLayer-1)%6].push_back(fdchit);
90  }
91  wire=fdchit->element;
92  layer=fdchit->gLayer;
93  }
94 
95  // Loop over packages, creating intersection points
96  for(unsigned int package=1; package<=4; package++){
97  MakeIntersectionPoints(fdchits_by_package[package-1]);
98  //MakeRestrictedIntersectionPoints(fdchits_by_package[package-1]);
99  }
100 
101  return NOERROR;
102 }
103 
104 //------------------
105 // MakeIntersectionPoints
106 //------------------
107 void DFDCIntersection_factory::MakeIntersectionPoints(vector<vector<const DFDCHit*> >&hits_by_layer)
108 {
109  // hits_by_layer should contain all of the wire hits in a single
110  // FDC package. They are stored in set of nested vectors with the
111  // outer vector being the list of layers (always 6 in length)
112  // and the inner vector the hits within the layer.
113 
114  // Loop over layers
115  for(unsigned int i=0; i<hits_by_layer.size()-1; i++){
116  vector<const DFDCHit*> &layer1 = hits_by_layer[i];
117  vector<const DFDCHit*> &layer2 = hits_by_layer[i+1];
118 
119  FindIntersections(layer1, layer2, _data);
120  }
121 }
122 
123 //------------------
124 // MakeRestrictedIntersectionPoints
125 //------------------
126 void DFDCIntersection_factory::MakeRestrictedIntersectionPoints(vector<vector<const DFDCHit*> >&hits_by_layer)
127 {
128  // hits_by_layer should contain all of the wire hits in a single
129  // FDC package. They are stored in a set of nested vectors with the
130  // outer vector being the list of layers (always 6 in length)
131  // and the inner vector the hits within the layer.
132 
133  // This method will first make a list of all the intersection points
134  // between adjacent layers in the package. It will then make a list
135  // which intersection points to keep by finding those that meet
136  // one of the following criteria:
137  //
138  // 1. For one of the wires used, this is the only intersection
139  // point at this z location.
140  //
141  // 2. The intersection point has a corresponding intersection
142  // point in an adjacent plane (i.e. a triple wire coincidence)
143  //
144  // Note that there is a class of hits this method fails to accept
145  // at the moment. This is when 2 wire planes have 2 hits each leading
146  // to 4 intersections, only 2 of which are real. If the adjacent
147  // plane has only one hit that is in coincidence with only
148  // 1 of the points, the other "true" point could be inferred. At this
149  // point, the second true point will not be kept in these cases
150  // (or similar ones with more tracks).
151 
152  // First, get all of the intersection points
153  vector<vector<DFDCIntersection *> > intersections(hits_by_layer.size()-1);
154  for(unsigned int i=0; i<hits_by_layer.size()-1; i++){
155 
156  // Find intersection points between the current layer and next one downstream
157  FindIntersections(hits_by_layer[i], hits_by_layer[i+1], intersections[i]);
158  }
159 
160  // Here we need to make 2 lists for every wire hit. The first is of
161  // intersection points made from intersections with its upstream
162  // neighbor and the second from intersections with its downstream
163  // neighbor.
164  map<const DFDCWire*, vector<DFDCIntersection*> > upstr;
165  map<const DFDCWire*, vector<DFDCIntersection*> > dnstr;
166  for(unsigned int i=0; i<intersections.size(); i++){
167  vector<DFDCIntersection *> &layer = intersections[i];
168 
169  for(unsigned int j=0; j<layer.size(); j++){
170  DFDCIntersection *fdcinter = layer[j];
171  // Note here that wire2 is in the downstream layer for the hit
172  // and wire1 is the upstream layer. Therefore, this intersection
173  // should be upstream for wire2 and downstream for wire 1.
174  upstr[fdcinter->wire2].push_back(fdcinter);
175  dnstr[fdcinter->wire1].push_back(fdcinter);
176  }
177  }
178 
179  // Keep track of the "good" intersection objects in a map.
180  map<DFDCIntersection *, bool> intersects_to_keep;
181 
182  // Loop over both the upstr and dnstr lists to find wires
183  // with only one hit in either its upstream or downstream
184  // plane.
185  map<const DFDCWire*, vector<DFDCIntersection*> >::iterator iter;
186  for(iter=upstr.begin(); iter!=upstr.end(); iter++){
187  vector<DFDCIntersection*> &hits = iter->second;
188  if(hits.size()==1)intersects_to_keep[hits[0]] = true;
189  }
190  for(iter=dnstr.begin(); iter!=dnstr.end(); iter++){
191  vector<DFDCIntersection*> &hits = iter->second;
192  if(hits.size()==1)intersects_to_keep[hits[0]] = true;
193  }
194 
195  // Find triple wire plane coincidences
196  // Loop over intersection layers
197  for(unsigned int i=0; i<intersections.size()-1; i++){
198 
199  // Loop over intersections in this intersection layer
200  vector<DFDCIntersection *> &layer1 = intersections[i];
201  for(unsigned int j=0; j<layer1.size(); j++){
202  DFDCIntersection *int1 = layer1[j];
203 
204  // Loop over hits in downstream layer
205  vector<DFDCIntersection *> &layer2 = intersections[i+1];
206  for(unsigned int k=0; k<layer2.size(); k++){
207  DFDCIntersection *int2 = layer2[k];
208  // Only interested in hits that share a wire
209  if(int1->wire2 != int2->wire1)continue;
210 
211  // Find distance between hits
212  double dx = int1->pos.X() - int2->pos.X();
213  double dy = int1->pos.Y() - int2->pos.Y();
214  double dist2 = dx*dx + dy*dy;
215 
216  if(dist2<MAX_DIST2){
217  intersects_to_keep[int1] = true;
218  intersects_to_keep[int2] = true;
219  }
220  }
221  }
222  }
223 
224  // Loop over all intersection points. Any that are in "intersects_to_keep"
225  // copy to _data. All others, delete.
226  for(unsigned int i=0; i<intersections.size(); i++){
227  for(unsigned int j=0; j<intersections[i].size(); j++){
228  DFDCIntersection *intersection = intersections[i][j];
229  bool keep = intersects_to_keep[intersection];
230  if(keep){
231  _data.push_back(intersection);
232  }else{
233  delete intersection;
234  }
235  }
236  }
237 }
238 
239 //------------------
240 // FindIntersections
241 //------------------
242 void DFDCIntersection_factory::FindIntersections(vector<const DFDCHit*> &layer1, vector<const DFDCHit*> &layer2, vector<DFDCIntersection*> &intersections)
243 {
244  // Loop over hits in layer1
245  for(unsigned int j=0; j<layer1.size(); j++){
246  const DFDCHit* hit1 = layer1[j];
247  const DFDCWire *wire1 =fdcwires[hit1->gLayer-1][hit1->element-1];
248  if(!wire1){
249  _DBG_<<"No wire for layer="<<hit1->gLayer<<" wire="<<hit1->element<<" !!"<<endl;
250  continue;
251  }
252  DVector2 a(wire1->origin.X(), wire1->origin.Y());
253  DVector2 b(wire1->udir.X(), wire1->udir.Y());
254  b /= b.Mod();
255 
256  // Loop over hits in layer2
257  for(unsigned int k=0; k<layer2.size(); k++){
258  const DFDCHit* hit2 = layer2[k];
259  const DFDCWire *wire2 =fdcwires[hit2->gLayer-1][hit2->element-1];
260  if(!wire2){
261  _DBG_<<"No wire for layer="<<hit2->gLayer<<" wire="<<hit2->element<<" !!"<<endl;
262  continue;
263  }
264  DVector2 c(wire2->origin.X(), wire2->origin.Y());
265  DVector2 d(wire2->udir.X(), wire2->udir.Y());
266  d /= d.Mod();
267 
268  // Find intersection of the projections of wires 1 and 2 on the X/Y plane
269  DVector2 e = a - c;
270  double gamma = d*b;
271  double beta = (d*e - gamma*(b*e))/(1.0-gamma*gamma);
272  if(!isfinite(beta))continue; // wires must be parallel
273  if(fabs(beta) > wire2->L/2.0)continue; // intersection is past end of wire
274  double alpha = beta*gamma - b*e;
275  if(fabs(alpha) > wire1->L/2.0)continue; // intersection is past end of wire
276 
277  DVector2 x = c + beta*d; // intersection point in X/Y
278 
279  // Add the intersection point
280  DFDCIntersection *fdcint = new DFDCIntersection;
281  fdcint->hit1 = hit1;
282  fdcint->hit2 = hit2;
283  fdcint->wire1 = wire1;
284  fdcint->wire2 = wire2;
285  fdcint->pos.SetXYZ(x.X(), x.Y(), (wire1->origin.Z()+wire2->origin.Z())/2.0);
286 
287  intersections.push_back(fdcint);
288  }
289  }
290 }
291 
DApplication * dapp
void MakeRestrictedIntersectionPoints(vector< vector< const DFDCHit * > > &hits_by_layer)
int type
Definition: DFDCHit.h:44
const DFDCWire * wire1
Int_t layer
TVector2 DVector2
Definition: DVector2.h:9
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
const DFDCHit * hit1
#define c
static vector< vector< DFDCWire * > > fdcwires
jerror_t brun(JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
int gLayer
Definition: DFDCHit.h:28
void MakeIntersectionPoints(vector< vector< const DFDCHit * > > &hits_by_layer)
jerror_t init(void)
Called once at program start.
const DFDCHit * hit2
TEllipse * e
const double alpha
DGeometry * GetDGeometry(unsigned int run_number)
#define _DBG_
Definition: HDEVIO.h:12
int element
Definition: DFDCHit.h:25
bool GetFDCWires(vector< vector< DFDCWire * > > &fdcwires) const
Definition: DGeometry.cc:1078
const DFDCWire * wire2
void FindIntersections(vector< const DFDCHit * > &layer1, vector< const DFDCHit * > &layer2, vector< DFDCIntersection * > &intersections)
class DFDCHit: definition for a basic FDC hit data type.
Definition: DFDCHit.h:20
jerror_t evnt(JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.