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);
41 _DBG_<<
"FDC geometry not available!" <<endl;
56 for(
int package=1; package<=4; package++){
58 fdchits_by_package[package-1][
layer-1].clear();
63 vector<const DFDCHit*> fdchits;
69 if(fdchits.size()>(5+5+1)*24*4){
71 _DBG_<<
"Too many hits in FDC! Intersection point reconstruction in FDC bypassed for event "<<loop->GetJEvent().GetEventNumber()<<endl;
78 for(
unsigned int i=0; i<fdchits.size(); i++){
79 const DFDCHit *fdchit = fdchits[i];
81 if(fdchit->
type != 0)
continue;
83 _DBG_<<
"FDC gLayer out of range! ("<<fdchit->
gLayer<<
" is not between 1 and 24)"<<endl;
88 int package = (fdchit->gLayer-1)/6 + 1;
89 fdchits_by_package[package-1][(fdchit->
gLayer-1)%6].push_back(fdchit);
96 for(
unsigned int package=1; package<=4; package++){
97 MakeIntersectionPoints(fdchits_by_package[package-1]);
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];
119 FindIntersections(layer1, layer2, _data);
153 vector<vector<DFDCIntersection *> > intersections(hits_by_layer.size()-1);
154 for(
unsigned int i=0; i<hits_by_layer.size()-1; i++){
157 FindIntersections(hits_by_layer[i], hits_by_layer[i+1], intersections[i]);
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];
169 for(
unsigned int j=0; j<layer.size(); j++){
174 upstr[fdcinter->
wire2].push_back(fdcinter);
175 dnstr[fdcinter->
wire1].push_back(fdcinter);
180 map<DFDCIntersection *, bool> intersects_to_keep;
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;
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;
197 for(
unsigned int i=0; i<intersections.size()-1; i++){
200 vector<DFDCIntersection *> &layer1 = intersections[i];
201 for(
unsigned int j=0; j<layer1.size(); j++){
205 vector<DFDCIntersection *> &layer2 = intersections[i+1];
206 for(
unsigned int k=0; k<layer2.size(); k++){
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;
217 intersects_to_keep[int1] =
true;
218 intersects_to_keep[int2] =
true;
226 for(
unsigned int i=0; i<intersections.size(); i++){
227 for(
unsigned int j=0; j<intersections[i].size(); j++){
229 bool keep = intersects_to_keep[intersection];
231 _data.push_back(intersection);
245 for(
unsigned int j=0; j<layer1.size(); j++){
246 const DFDCHit* hit1 = layer1[j];
257 for(
unsigned int k=0; k<layer2.size(); k++){
258 const DFDCHit* hit2 = layer2[k];
271 double beta = (d*e - gamma*(b*
e))/(1.0-gamma*gamma);
272 if(!isfinite(beta))
continue;
273 if(fabs(beta) > wire2->
L/2.0)
continue;
274 double alpha = beta*gamma - b*
e;
275 if(fabs(alpha) > wire1->
L/2.0)
continue;
283 fdcint->
wire1 = wire1;
284 fdcint->
wire2 = wire2;
285 fdcint->
pos.SetXYZ(x.X(), x.Y(), (wire1->
origin.Z()+wire2->
origin.Z())/2.0);
287 intersections.push_back(fdcint);
void MakeRestrictedIntersectionPoints(vector< vector< const DFDCHit * > > &hits_by_layer)
static vector< vector< DFDCWire * > > fdcwires
jerror_t brun(JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
void MakeIntersectionPoints(vector< vector< const DFDCHit * > > &hits_by_layer)
jerror_t init(void)
Called once at program start.
DGeometry * GetDGeometry(unsigned int run_number)
bool GetFDCWires(vector< vector< DFDCWire * > > &fdcwires) const
void FindIntersections(vector< const DFDCHit * > &layer1, vector< const DFDCHit * > &layer2, vector< DFDCIntersection * > &intersections)
class DFDCHit: definition for a basic FDC hit data type.
jerror_t evnt(JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.