Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DTrackCandidate_factory_FDC.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: DTrackCandidate_factory.cc
4 // Created: Mon Jul 18 15:23:04 EDT 2005
5 // Creator: davidl (on Darwin wire129.jlab.org 7.8.0 powerpc)
6 //
7 
8 #include <iostream>
9 #include <iomanip>
10 using namespace std;
11 
12 #include <TROOT.h>
13 
14 #include <math.h>
15 
17 #include "DANA/DApplication.h"
18 #include "DQuickFit.h"
19 #include "JANA/JGeometry.h"
21 #include "DVector2.h"
22 #include "FDC/DFDCGeometry.h"
23 #include "DHoughFind.h"
24 
25 
27  return hit1->hit->pos.Z() < hit2->hit->pos.Z();
28 }
29 
30 
31 //------------------
32 // DTrackCandidate_factory_FDC
33 //------------------
35 {
36 #if 0
37  // Set defaults
38  MAX_SEED_DIST = 5.0;
39  gPARMS->SetDefaultParameter("TRKFIND:MAX_SEED_DIST", MAX_SEED_DIST);
40 #endif
41 }
42 
43 //------------------
44 // init
45 //------------------
47 {
48  TARGET_Z_MIN = 65.0 - 15.0;
49  TARGET_Z_MAX = 65.0 + 15.0;
50 
51  MAX_HIT_DIST = 5.0;
52  MAX_HIT_DIST2 = MAX_HIT_DIST*MAX_HIT_DIST;
53 
54  return NOERROR;
55 }
56 
57 //------------------
58 // brun
59 //------------------
60 jerror_t DTrackCandidate_factory_FDC::brun(JEventLoop *loop, int32_t runnumber)
61 {
62  return NOERROR;
63 }
64 
65 //------------------
66 // fini
67 //------------------
69 {
70  return NOERROR;
71 }
72 
73 //------------------
74 // evnt
75 //------------------
76 jerror_t DTrackCandidate_factory_FDC::evnt(JEventLoop *loop, uint64_t eventnumber)
77 {
78 
79  // Get the hits into the trkhits vector
80  GetTrkHits(loop);
81 
82  // Find seeds from X/Y projections
83  vector<DFDCSeed> seeds;
84  FindSeeds(seeds);
85 
86  // Loop over seeds and fit in phi-z plane to find z and theta
87  for(unsigned int i=0; i<seeds.size(); i++){
88  DFDCSeed &seed = seeds[i];
89  if(debug_level>3)_DBG_<<"----- Seed "<<i<<" ------"<<endl;
90  if(!seed.valid)continue;
91 
92  // Fit seed hits to get theta and vertex z position
93  //FindThetaZ(seed);
94  //if(!seed.valid)continue;
95 
96  // copy fit results to local variables (makes it easier to debug)
97  double p_trans = seed.p_trans;
98  double phi = seed.phi;
99  double q = seed.q;
100  double theta = seed.theta;
101  double z_vertex = seed.z_vertex;
102 
103  if(debug_level>3)_DBG_<<"p_trans="<<p_trans<<" phi="<<phi<<" theta="<<theta<<" p="<<p_trans/sin(theta)<<endl;
104 
105  //Make a track candidate from results
106  DTrackCandidate *can = new DTrackCandidate;
107  DVector3 pos, mom;
108  pos.SetXYZ(0.0, 0.0, z_vertex);
109  mom.SetMagThetaPhi(p_trans/sin(theta), theta, phi);
110  can->setPosition(pos);
111  can->setMomentum(mom);
112  can->setPID((q > 0.0) ? PiPlus : PiMinus);
113 
114  _data.push_back(can);
115  }
116 
117  return NOERROR;
118 }
119 
120 //------------------
121 // GetTrkHits
122 //------------------
124 {
125  // Clear out old hits
126  for(unsigned int i=0; i<fdctrkhits.size(); i++){
127  delete fdctrkhits[i];
128  }
129  fdctrkhits.clear();
130 
131  // Get hits
132  vector<const DFDCIntersection*> fdcintersects;
133  loop->Get(fdcintersects);
134 
135  // Create a DFDCTrkHit object for each DFDCIntersection object
136  for(unsigned int i=0; i<fdcintersects.size(); i++){
137 
138  DFDCTrkHit *trkhit = new DFDCTrkHit;
139  trkhit->hit = fdcintersects[i];
140  trkhit->flags = NOISE;
141 
142  fdctrkhits.push_back(trkhit);
143  }
144 
145  // Filter out noise hits. All hits are initially flagged as "noise".
146  // Hits with a neighbor within MAX_HIT_DIST have their noise flags cleared.
147  for(unsigned int i=0; i<fdctrkhits.size(); i++){ // cut above should ensure cdctrkhits.size() is at least 1
148  DFDCTrkHit *trkhit1 = fdctrkhits[i];
149  for(unsigned int j=i+1; j<fdctrkhits.size(); j++){
150  DFDCTrkHit *trkhit2 = fdctrkhits[j];
151  if(!(trkhit2->flags & NOISE))continue; // this hit already has noise flag cleared
152  double d2 = trkhit1->Dist2(trkhit2);
153  double deltaz = fabs(trkhit1->hit->pos.Z() - trkhit2->hit->pos.Z());
154  if(debug_level>4)_DBG_<<" -- Dist hits "<<i<<" and "<<j<<" : deltaR="<<sqrt(d2)<<" deltaZ="<<deltaz<<endl;
155  if(d2<MAX_HIT_DIST2 && deltaz<12.0){
156  trkhit1->flags &= ~NOISE;
157  trkhit2->flags &= ~NOISE;
158  }// if
159  }// j
160  }// i
161 }
162 
163 
164 
165 //------------------
166 // FindSeeds
167 //------------------
168 void DTrackCandidate_factory_FDC::FindSeeds(vector<DFDCSeed> &seeds)
169 {
170  // Create a DHoughFind object to find a seed via the Hough Transform method
171  // We create it once here to avoid the overhead of reallocating memory
172  // each time through the loop below.
173  //DHoughFind hough(-400.0, +400.0, -400.0, +400.0, 100, 100);
174 
175  // Loop until we find all seeds
176  unsigned int last_available_hits = 0;
177  for(int i=0; i<12; i++){ // limit the number of iterations in this loop
178  // Check if we should exit the loop due to lack of available
179  // hits. This is the primary point of exit for this loop.
180  unsigned int Navailable_hits = NumAvailableHits();
181  if(debug_level>3)_DBG_<<"Navailable_hits="<<Navailable_hits<<" last_available_hits="<<last_available_hits<<endl;
182  if(Navailable_hits<3)break;
183  if(Navailable_hits==last_available_hits)break;
184  last_available_hits = Navailable_hits;
185 
186  // Set the limits for the rough Hough histogram
187  hough.SetLimits(-400.0, +400.0, -400.0, +400.0, 100, 100);
188  hough.ClearPoints();
189 
190  // Add all points not already marked as unavailable
191  for(unsigned int i=0; i<fdctrkhits.size(); i++){
192  DFDCTrkHit *fdctrkhit = fdctrkhits[i];
193  if(fdctrkhit->flags&USED)continue;
194  if(fdctrkhit->flags&NOISE)continue;
195  if(fdctrkhit->flags&OUT_OF_TIME)continue;
196  const DFDCIntersection *hit = fdctrkhit->hit;
197 
198  hough.AddPoint(hit->pos.X(), hit->pos.Y());
199  }
200 
201  DVector2 Ro = hough.Find();
202  if(debug_level>3)_DBG_<<"Rough: GetMaxBinContent="<<hough.GetMaxBinContent()<<" x0="<<Ro.X()<<" y0="<<Ro.Y()<<endl;
203  if(debug_level>6)hough.PrintHist();
204  if(hough.GetMaxBinContent()<10.0)break;
205 
206  // Zoom in on resonanace a little
207  double width = 60.0;
208  hough.SetLimits(Ro.X()-width, Ro.X()+width, Ro.Y()-width, Ro.Y()+width, 100, 100);
209  Ro = hough.Find();
210  if(debug_level>3)_DBG_<<"Medium: GetMaxBinContent="<<hough.GetMaxBinContent()<<" x0="<<Ro.X()<<" y0="<<Ro.Y()<<endl;
211  if(debug_level>5)hough.PrintHist();
212 
213  // Zoom in on resonanace once more
214  width = 8.0;
215  hough.SetLimits(Ro.X()-width, Ro.X()+width, Ro.Y()-width, Ro.Y()+width, 100, 100);
216  Ro = hough.Find();
217  if(debug_level>3)_DBG_<<"Fine: GetMaxBinContent="<<hough.GetMaxBinContent()<<" x0="<<Ro.X()<<" y0="<<Ro.Y()<<endl;
218  if(debug_level>5)hough.PrintHist();
219 
220  DFDCSeed seed;
221 
222  seed.valid = true;
223  seed.q = +1.0;
224  seed.z_vertex = 65.0;
225  seed.theta = M_PI/2.0;
226 
227  seed.r0 = Ro.Mod();
228  seed.x0 = Ro.X();
229  seed.y0 = Ro.Y();
230  seed.p_trans = 0.003*2.1*seed.r0;
231 
232  // Set the VALID_HIT flag on hits consistent with this circle
233  // and calculate the phi angle relative to the center of the
234  // circle for the valid hits.
235  FillSeedHits(seed);
236  if(seed.hits.size()<3)continue; // require at least 5 hits
237 
238  // Set phi here since FillSeedHits may flip the sign of seed.q
239  seed.phi = Ro.Phi() - seed.q*M_PI/2.0;
240 
241  // Find the theta and z-vertex of the seed and flag any used hits as USED
242  FindThetaZ(seed);
243 
244  // Add to list of seeds (i.e. tracks)
245  seeds.push_back(seed);
246  }
247 }
248 
249 //------------------
250 // FillSeedHits
251 //------------------
253 {
254  // Loop over available hits and add any consistent with the circle
255  // parameters to the list of hits for this seed.
256  DVector2 Ro(seed.x0, seed.y0);
257  seed.hits.clear();
258  for(unsigned int i=0; i<fdctrkhits.size(); i++){
259  DFDCTrkHit *fdctrkhit = fdctrkhits[i];
260  fdctrkhit->flags &= ~ON_CIRCLE; // clear ON_CIRCLE flag, it gets set below if applicable
261  if(fdctrkhit->flags&USED)continue;
262  if(fdctrkhit->flags&NOISE)continue;
263  if(fdctrkhit->flags&OUT_OF_TIME)continue;
264  const DFDCIntersection *hit = fdctrkhit->hit;
265 
266  // Calculate distance between Hough transformed line (i.e.
267  // the line on which a circle that passes through both the
268  // origin and the point at hit->pos) and the circle center.
269  DVector2 h(hit->pos.X()/2.0, hit->pos.Y()/2.0);
270  DVector2 g(h.Y(), -h.X());
271  g /= g.Mod();
272  DVector2 Ro_minus_h(seed.x0-h.X(), seed.y0-h.Y());
273  //Ro_minus_h /= Ro_minus_h.Mod();
274  double dist = fabs(g.X()*Ro_minus_h.Y() - g.Y()*Ro_minus_h.X());
275 
276  // If this is not close enough to the found circle's center,
277  // reject it for this seed.
278  if(debug_level>3)_DBG_<<"dist="<<dist<<endl;
279  if(dist > 2.0){
280  fdctrkhit->flags &= ~VALID_HIT;
281  continue;
282  }
283 
284  fdctrkhit->flags |= ON_CIRCLE | VALID_HIT;
285  seed.hits.push_back(fdctrkhit);
286  }
287 
288  // Sort hits by increasing z
289  stable_sort(seed.hits.begin(), seed.hits.end(), FDCSortByZincreasing);
290 
291  // Below, we fill in the phi_hit field of the DFDCTrkHit objects for
292  // this candidate. We do so incrementally to try and accomodate
293  // phi values greater than 2pi. Initialize the "last" direction as
294  // pointing back to the beamline.
295  DVector2 last_dir = -1.0*Ro/Ro.Mod();
296  double last_phi = 0.0;
297 
298  // Loop over hits, filling in phi_hit
299  int Nq = 0;
300  for(unsigned int i=0; i<seed.hits.size(); i++){
301  DFDCTrkHit *fdctrkhit = seed.hits[i];
302  const DVector3 &pos = fdctrkhit->hit->pos;
303 
304  // Calculate phi. We do this trivially for now
305  DVector2 p(pos.X() , pos.Y());
306  DVector2 p_minus_Ro = p - Ro;
307  p_minus_Ro/=p_minus_Ro.Mod();
308  double delta_phi = p_minus_Ro.Phi() - last_dir.Phi();
309  while(delta_phi>+M_PI)delta_phi-=2.0*M_PI;
310  while(delta_phi<-M_PI)delta_phi+=2.0*M_PI;
311  fdctrkhit->phi_hit = last_phi + delta_phi;
312  last_phi = fdctrkhit->phi_hit;
313  last_dir = p_minus_Ro;
314  if(debug_level>3)_DBG_<<"phi_hit="<<fdctrkhit->phi_hit<<" z="<<fdctrkhit->hit->pos.Z()<<" phi_hit/z="<<fdctrkhit->phi_hit/(fdctrkhit->hit->pos.Z()-65.0)<<" delta_phi="<<delta_phi<<endl;
315  if(fdctrkhit->hit->wire2->layer<=6){
316  Nq += fdctrkhit->phi_hit>0.0 ? +1:-1;
317  }
318  }
319 
320  // If Nq has too few entries, then look to the second package
321  if(abs(Nq)<2){
322  for(unsigned int i=0; i<seed.hits.size(); i++){
323  DFDCTrkHit *fdctrkhit = seed.hits[i];
324  if(fdctrkhit->hit->wire2->layer<=12){
325  Nq += fdctrkhit->phi_hit>0.0 ? +1:-1;
326  }
327  }
328  }
329 
330  if(Nq<0)seed.q = -seed.q;
331 }
332 
333 //------------------
334 // NumAvailableHits
335 //------------------
337 {
338  // Loop over all hits and count the ones that are still
339  // available for making a new seed.
340  unsigned int N=0;
341  for(unsigned int i=0; i<fdctrkhits.size(); i++){
342  DFDCTrkHit *fdctrkhit = fdctrkhits[i];
343  if(fdctrkhit->flags&USED)continue;
344  if(fdctrkhit->flags&NOISE)continue;
345  if(fdctrkhit->flags&OUT_OF_TIME)continue;
346  if(fdctrkhit->flags&CANT_BE_IN_SEED)continue;
347  N++;
348  }
349 
350  return N;
351 }
352 
353 //------------------
354 // FindThetaZ
355 //------------------
357 {
358  FindTheta(seed, TARGET_Z_MIN, TARGET_Z_MAX);
359  FindZ(seed, seed.theta_min, seed.theta_max);
360 
361  // If z_vertex is not inside the target limits, then flag this
362  // seed as invalid.
363  if(seed.z_vertex<TARGET_Z_MIN || seed.z_vertex>TARGET_Z_MAX){
364  if(debug_level>3)_DBG_<<"Seed z-vertex outside of target range (z="<<seed.z_vertex<<" TARGET_Z_MIN="<<TARGET_Z_MIN<<" TARGET_Z_MAX="<<TARGET_Z_MAX<<endl;
365  seed.valid=false;
366  }
367 
368  // If the seed is valid, mark all hits that are consistent with the track
369  // as USED
370  if(seed.valid){
371  for(unsigned int i=0; i<seed.hits.size(); i++){
372  DFDCTrkHit *trkhit = seed.hits[i];
373  if(trkhit->flags&IN_Z_RANGE)trkhit->flags |= USED;
374  }
375  }else{
376  // If the seed is not valid, then we need to mark at least one of the
377  // hits us unusable on the next seed. Otherwise, the same seed will be found
378  // over and over and over... We just mark the first hit in the seed.
379  for(unsigned int i=0; i<seed.hits.size(); i++){
380  DFDCTrkHit *trkhit = seed.hits[i];
381  if(trkhit->flags&ON_CIRCLE)trkhit->flags |= CANT_BE_IN_SEED;
382  }
383  }
384 
385  return;
386 }
387 
388 //------------------
389 // FindTheta
390 //------------------
391 void DTrackCandidate_factory_FDC::FindTheta(DFDCSeed &seed, double target_z_min, double target_z_max)
392 {
393  /// Find the theta value using the hits from <i>seed</i>.
394  /// The value of seed.r0 is used to calculate theta.
395  ///
396  /// This uses a histogramming technique that looks at the overlaps of the
397  /// angle ranges subtended by each hit between the given target limits.
398  /// The overlaps usually lead to a range of values for theta. The limits
399  /// of these are stored in the theta_min and theta_max fields of the seed.
400  /// The centroid of the range is stored in the theta field.
401 
402  // We use a simple array to store our histogram here. We don't want to use
403  // ROOT histograms because they are not thread safe.
404  unsigned int Nbins = 1000;
405  unsigned int hist[Nbins];
406  for(unsigned int i=0; i<Nbins; i++)hist[i] = 0; // clear histogram
407  double bin_width = 2.0*M_PI/(double)Nbins;
408  double hist_low_limit = -M_PI; // lower edge of histogram limits
409 
410  // Loop over valid hits, filling the histogram
411  double &r0 = seed.r0;
412  for(unsigned int i=0; i<seed.hits.size(); i++){
413  DFDCTrkHit *trkhit = seed.hits[i];
414  if(!trkhit->flags&VALID_HIT)continue;
415  if(!trkhit->flags&ON_CIRCLE)continue;
416 
417  // Calculate upper and lower limits in theta
418  double alpha = r0*trkhit->phi_hit;
419  if(seed.q<0.0)alpha = -alpha;
420  double z_hit = trkhit->hit->pos.Z();
421  double tmin = atan2(alpha, z_hit - target_z_min);
422  double tmax = atan2(alpha, z_hit - target_z_max);
423  if(tmin>tmax){
424  double tmp = tmin;
425  tmin=tmax;
426  tmax=tmp;
427  }
428  if(debug_level>3)_DBG_<<" -- phi_hit="<<trkhit->phi_hit<<" z_hit="<<z_hit<<endl;
429  if(debug_level>3)_DBG_<<" -- tmin="<<tmin<<" tmax="<<tmax<<endl;
430 
431  // Copy theta limits into trkhit
432  trkhit->theta_min = tmin;
433  trkhit->theta_max = tmax;
434 
435  // Find index of bins corresponding to tmin and tmax
436  unsigned int imin = (unsigned int)floor((tmin-hist_low_limit)/bin_width);
437  unsigned int imax = (unsigned int)floor((tmax-hist_low_limit)/bin_width);
438 
439  // If entire range of this hit is outside of the histogram limit
440  // then discard this hit.
441  if(imin>=Nbins)continue;
442 
443  // Clip limits of bin range to our histogram limits
444  if(imin>=Nbins)imin=Nbins-1;
445  if(imax>=Nbins)imax=Nbins-1;
446 
447  // Increment all bins between imin and imax
448  for(unsigned int j=imin; j<=imax; j++)hist[j]++;
449  }
450 
451  // Look for the indexes of the plateau
452  unsigned int istart=0;
453  unsigned int iend=0;
454  for(unsigned int i=1; i<Nbins; i++){
455  if(hist[i]>hist[istart]){
456  istart = i;
457  if(debug_level>3)_DBG_<<" -- istart="<<istart<<" (theta="<<hist_low_limit + bin_width*(0.5+(double)istart)<<" , N="<<hist[i]<<")"<<endl;
458  }
459  if(hist[i] == hist[istart])iend = i;
460  }
461 
462  // If there are no entries in the histogram, then flag this seed as invalid
463  if(hist[istart]==0.0){
464  if(debug_level>3)_DBG_<<" - No entries in theta hist. Seed aborted."<<endl;
465  seed.valid=false;
466  }
467 
468  // Calculate theta limits
469  seed.theta_min = hist_low_limit + bin_width*(0.5+(double)istart);
470  seed.theta_max = hist_low_limit + bin_width*(0.5+(double)iend);
471  seed.theta = (seed.theta_max + seed.theta_min)/2.0;
472  if(debug_level>3)_DBG_<<"istart="<<istart<<" iend="<<iend<<" theta_min="<<seed.theta_min<<" theta_max="<<seed.theta_max<<endl;
473 
474  // Mark all of the on circle hits that have a theta range consistent with seed.theta
475  // as being IN_THETA_RANGE.
476  for(unsigned int i=0; i<seed.hits.size(); i++){
477  DFDCTrkHit *trkhit = seed.hits[i];
478  trkhit->flags &= ~IN_THETA_RANGE;
479  if(!trkhit->flags&VALID_HIT)continue;
480  if(!trkhit->flags&ON_CIRCLE)continue;
481  if(trkhit->theta_min > seed.theta)continue;
482  if(trkhit->theta_max < seed.theta)continue;
483  trkhit->flags |= IN_THETA_RANGE;
484  }
485 
486  // If theta is negative, then we probably chose the wrong sign. Flip
487  // it now.
488  if(seed.theta<0.0){
489  seed.theta_min *= -1.0;
490  seed.theta_max *= -1.0;
491  seed.theta *= -1.0;
492  seed.q *= -1.0;
493  seed.phi += M_PI;
494  if(seed.phi > 2.0*M_PI)seed.phi -= 2.0*M_PI;
495  }
496 }
497 
498 //------------------
499 // FindZ
500 //------------------
501 void DTrackCandidate_factory_FDC::FindZ(DFDCSeed &seed, double theta_min, double theta_max)
502 {
503  /// Find the z value of the vertex using the valid stereo hits from <i>seed</i>. The values
504  /// for phi_hit and pos.Z() are assumed to be valid as is the status of the
505  /// VALID_HIT bit in flags.
506  ///
507  /// This uses a histogramming technique that looks at the overlaps of the
508  /// z ranges subtended by each hit between the given theta limits.
509  /// The overlaps usually lead to a range of values for z_vertex. The limits
510  /// of these are stored in the z_min and z_max fields of the seed.
511  /// The centroid of the range is stored in the z_vertex field.
512 
513  // We use a simple array to store our histogram here. We don't want to use
514  // ROOT histograms because they are not thread safe.
515  unsigned int Nbins = 300;
516  unsigned int hist[Nbins];
517  for(unsigned int i=0; i<Nbins; i++)hist[i] = 0; // clear histogram
518  double bin_width = 0.5; // bins are 5mm
519  double hist_low_limit = 0.0; // lower edge of histogram limits
520 
521  // We effectively extend the theta_min and theta_max angles here
522  // a bit to include some error. The motivation is that if
523  // theta_min == theta_max that leads to z_min == z_max so there
524  // is little or no overlap of the z ranges of separate hits.
525  // For now, we hardwire this to 1 degree
526  double theta_err = 1.0/57.3;
527 
528  // Loop over valid hits, filling the histogram
529  double r0 = seed.r0;
530  double tan_alpha_max = tan(theta_max + theta_err)/r0;
531  double tan_alpha_min = tan(theta_min - theta_err)/r0;
532  if(tan_alpha_min<0.0)tan_alpha_min=0.0;
533  for(unsigned int i=0; i<seed.hits.size(); i++){
534  DFDCTrkHit *trkhit = seed.hits[i];
535  if(!trkhit->flags&VALID_HIT)continue;
536  if(!trkhit->flags&ON_CIRCLE)continue;
537  if(!trkhit->flags&IN_THETA_RANGE)continue;
538 
539  // Calculate upper and lower limits in z
540  double q_sign = seed.q>0.0 ? +1.0:-1.0;
541  double z_hit = trkhit->hit->pos.Z();
542  double zmin = z_hit - q_sign*trkhit->phi_hit/tan_alpha_min;
543  double zmax = z_hit - q_sign*trkhit->phi_hit/tan_alpha_max;
544  if(zmin>zmax){
545  double tmp = zmin;
546  zmin=zmax;
547  zmax=tmp;
548  }
549 
550  // Copy z limits into trkhit
551  trkhit->zmin = zmin;
552  trkhit->zmax = zmax;
553 
554  // Find index of bins corresponding to tmin and tmax
555  unsigned int imin = zmin<hist_low_limit ? 0:(unsigned int)floor((zmin-hist_low_limit)/bin_width);
556  unsigned int imax = zmax<hist_low_limit ? 0:(unsigned int)floor((zmax-hist_low_limit)/bin_width);
557 
558  if(debug_level>3)_DBG_<<" -- phi_hit="<<trkhit->phi_hit<<" z_hit="<<z_hit<<endl;
559  if(debug_level>3)_DBG_<<" -- zmin="<<zmin<<" zmax="<<zmax<<endl;
560  if(debug_level>3)_DBG_<<" -- imin="<<imin<<" imax="<<imax<<endl;
561 
562  // If entire range of this hit is outside of the histogram limit
563  // then discard this hit.
564  if(imax<=0 || imin>=Nbins)continue;
565 
566  // Clip limits of bin range to our histogram limits
567  if(imin>=Nbins)imin=Nbins-1;
568  if(imax>=Nbins)imax=Nbins-1;
569 
570  // Increment all bins between imin and imax
571  for(unsigned int j=imin; j<=imax; j++)hist[j]++;
572  }
573 
574  // Look for the indexes of the plateau
575  unsigned int istart=(unsigned int)((TARGET_Z_MIN-hist_low_limit)/bin_width);
576  unsigned int iend=0;
577  for(unsigned int i=1; i<Nbins; i++){
578 
579  // Only look in Target area
580  double z = hist_low_limit + bin_width*(0.5+(double)i);
581  if(z<TARGET_Z_MIN || z>TARGET_Z_MAX)continue;
582 
583  if(hist[i]>hist[istart]){
584  istart = i;
585  if(debug_level>3)_DBG_<<" -- istart="<<istart<<" (z="<<hist_low_limit + bin_width*(0.5+(double)istart)<<" , N="<<hist[i]<<")"<<endl;
586  }
587  if(hist[i] == hist[istart])iend = i;
588  }
589 
590  // If there are no entries in the histogram, then flag this seed as invalid
591  if(hist[istart]==0.0){
592  if(debug_level>3)_DBG_<<" - No entries in z-vertex hist. Seed aborted."<<endl;
593  seed.valid=false;
594  }
595 
596  // Calculate z limits
597  seed.z_min = hist_low_limit + bin_width*(0.5+(double)istart);
598  seed.z_max = hist_low_limit + bin_width*(0.5+(double)iend);
599  seed.z_vertex = (seed.z_max + seed.z_min)/2.0;
600  if(debug_level>3)_DBG_<<"istart="<<istart<<" iend="<<iend<<" z_min="<<seed.z_min<<" z_max="<<seed.z_max<<" hits[istart]="<<hist[istart]<<endl;
601 
602  // Mark all of the hits that have a theta range consistent with seed.theta
603  // and a z_vertex consistent with seed.z_vertex as being IN_Z_RANGE.
604  for(unsigned int i=0; i<seed.hits.size(); i++){
605  DFDCTrkHit *trkhit = seed.hits[i];
606  trkhit->flags &= ~IN_Z_RANGE;
607  if(!trkhit->flags&VALID_HIT)continue;
608  if(!trkhit->flags&ON_CIRCLE)continue;
609  if(!trkhit->flags&IN_THETA_RANGE)continue;
610  if(trkhit->zmin > seed.z_vertex)continue;
611  if(trkhit->zmax < seed.z_vertex)continue;
612  trkhit->flags |= IN_Z_RANGE;
613  }
614 }
615 
void setMomentum(const DVector3 &aMomentum)
TVector2 DVector2
Definition: DVector2.h:9
TVector3 DVector3
Definition: DVector3.h:14
double zmax
virtual jerror_t brun(JEventLoop *loop, int32_t runnumber)
int layer
1-24
Definition: DFDCWire.h:16
void FindTheta(DFDCSeed &seed, double target_z_min, double target_z_max)
void FindSeeds(vector< DFDCSeed > &seeds)
const double alpha
void FindZ(DFDCSeed &seed, double theta_min, double theta_max)
#define _DBG_
Definition: HDEVIO.h:12
bool FDCSortByZincreasing(DTrackCandidate_factory_FDC::DFDCTrkHit *const &hit1, DTrackCandidate_factory_FDC::DFDCTrkHit *const &hit2)
virtual jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
Invoked via JEventProcessor virtual method.
virtual jerror_t fini(void)
Invoked via JEventProcessor virtual method.
&lt;A href=&quot;index.html#legend&quot;&gt; &lt;IMG src=&quot;CORE.png&quot; width=&quot;100&quot;&gt; &lt;/A&gt;
void setPID(Particle_t locPID)
double sqrt(double)
double sin(double)
void setPosition(const DVector3 &aPosition)
const DFDCWire * wire2
double zmin
int seed
TH1F * hist[Idx+1]
Definition: readhist.C:10