Bug Summary

File:libraries/TRACKING/DTrackCandidate_factory_FDC.cc
Location:line 148, column 25
Description:Dereference of null pointer

Annotated Source Code

1// $Id: DTrackCandidate_factory_FDC.cc 5619 2009-10-21 11:57:26Z staylor $
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>
10using namespace std;
11
12#include <TROOT.h>
13
14#include <math.h>
15
16#include "DTrackCandidate_factory_FDC.h"
17#include "DANA/DApplication.h"
18#include "DQuickFit.h"
19#include "JANA/JGeometry.h"
20#include "HDGEOMETRY/DMagneticFieldMap.h"
21#include "DVector2.h"
22#include "FDC/DFDCGeometry.h"
23#include "DHoughFind.h"
24
25
26bool FDCSortByZincreasing(DTrackCandidate_factory_FDC::DFDCTrkHit* const &hit1, DTrackCandidate_factory_FDC::DFDCTrkHit* const &hit2) {
27 return hit1->hit->pos.Z() < hit2->hit->pos.Z();
28}
29
30
31//------------------
32// DTrackCandidate_factory_FDC
33//------------------
34DTrackCandidate_factory_FDC::DTrackCandidate_factory_FDC()
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//------------------
46jerror_t DTrackCandidate_factory_FDC::init(void)
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//------------------
60jerror_t DTrackCandidate_factory_FDC::brun(JEventLoop *loop, int runnumber)
61{
62 return NOERROR;
63}
64
65//------------------
66// fini
67//------------------
68jerror_t DTrackCandidate_factory_FDC::fini(void)
69{
70 return NOERROR;
71}
72
73//------------------
74// evnt
75//------------------
76jerror_t DTrackCandidate_factory_FDC::evnt(JEventLoop *loop, int 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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
89<<" "<<"----- 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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
103<<" "<<"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->setCharge(q);
113
114 _data.push_back(can);
115 }
116
117 return NOERROR;
118}
119
120//------------------
121// GetTrkHits
122//------------------
123void DTrackCandidate_factory_FDC::GetTrkHits(JEventLoop *loop)
124{
125 // Clear out old hits
126 for(unsigned int i=0; i<fdctrkhits.size(); i++){
1
Loop condition is false. Execution continues on line 129
→
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++){
2
←
Loop condition is true. Entering loop body
→
3
←
Loop condition is false. Execution continues on line 147
→
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
4
←
Loop condition is true. Entering loop body
→
148 DFDCTrkHit *trkhit1 = fdctrkhits[i];
5
←
Dereference of null pointer
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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
154<<" "<<" -- 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//------------------
168void 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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
181<<" "<<"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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
202<<" "<<"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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
210<<" "<<"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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
217<<" "<<"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_PI3.14159265358979323846/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_PI3.14159265358979323846/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//------------------
252void DTrackCandidate_factory_FDC::FillSeedHits(DFDCSeed &seed)
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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
278<<" "<<"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 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_PI3.14159265358979323846)delta_phi-=2.0*M_PI3.14159265358979323846;
310 while(delta_phi<-M_PI3.14159265358979323846)delta_phi+=2.0*M_PI3.14159265358979323846;
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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
314<<" "<<"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//------------------
336unsigned int DTrackCandidate_factory_FDC::NumAvailableHits(void)
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//------------------
356void DTrackCandidate_factory_FDC::FindThetaZ(DFDCSeed &seed)
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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
364<<" "<<"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//------------------
391void 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_PI3.14159265358979323846/(double)Nbins;
408 double hist_low_limit = -M_PI3.14159265358979323846; // 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_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
428<<" "<<" -- phi_hit="<<trkhit->phi_hit<<" z_hit="<<z_hit<<endl;
429 if(debug_level>3)_DBG_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
429<<" "<<" -- 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(imax<0 || imin>=Nbins)continue;
442
443 // Clip limits of bin range to our histogram limits
444 if(imin<0)imin=0;
445 if(imin>=Nbins)imin=Nbins-1;
446 if(imax<0)imax=0;
447 if(imax>=Nbins)imax=Nbins-1;
448
449 // Increment all bins between imin and imax
450 for(unsigned int j=imin; j<=imax; j++)hist[j]++;
451 }
452
453 // Look for the indexes of the plateau
454 unsigned int istart=0;
455 unsigned int iend=0;
456 for(unsigned int i=1; i<Nbins; i++){
457 if(hist[i]>hist[istart]){
458 istart = i;
459 if(debug_level>3)_DBG_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
459<<" "<<" -- istart="<<istart<<" (theta="<<hist_low_limit + bin_width*(0.5+(double)istart)<<" , N="<<hist[i]<<")"<<endl;
460 }
461 if(hist[i] == hist[istart])iend = i;
462 }
463
464 // If there are no entries in the histogram, then flag this seed as invalid
465 if(hist[istart]==0.0){
466 if(debug_level>3)_DBG_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
466<<" "<<" - No entries in theta hist. Seed aborted."<<endl;
467 seed.valid=false;
468 }
469
470 // Calculate theta limits
471 seed.theta_min = hist_low_limit + bin_width*(0.5+(double)istart);
472 seed.theta_max = hist_low_limit + bin_width*(0.5+(double)iend);
473 seed.theta = (seed.theta_max + seed.theta_min)/2.0;
474 if(debug_level>3)_DBG_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
474<<" "<<"istart="<<istart<<" iend="<<iend<<" theta_min="<<seed.theta_min<<" theta_max="<<seed.theta_max<<endl;
475
476 // Mark all of the on circle hits that have a theta range consistent with seed.theta
477 // as being IN_THETA_RANGE.
478 for(unsigned int i=0; i<seed.hits.size(); i++){
479 DFDCTrkHit *trkhit = seed.hits[i];
480 trkhit->flags &= ~IN_THETA_RANGE;
481 if(!trkhit->flags&VALID_HIT)continue;
482 if(!trkhit->flags&ON_CIRCLE)continue;
483 if(trkhit->theta_min > seed.theta)continue;
484 if(trkhit->theta_max < seed.theta)continue;
485 trkhit->flags |= IN_THETA_RANGE;
486 }
487
488 // If theta is negative, then we probably chose the wrong sign. Flip
489 // it now.
490 if(seed.theta<0.0){
491 seed.theta_min *= -1.0;
492 seed.theta_max *= -1.0;
493 seed.theta *= -1.0;
494 seed.q *= -1.0;
495 seed.phi += M_PI3.14159265358979323846;
496 if(seed.phi > 2.0*M_PI3.14159265358979323846)seed.phi -= 2.0*M_PI3.14159265358979323846;
497 }
498}
499
500//------------------
501// FindZ
502//------------------
503void DTrackCandidate_factory_FDC::FindZ(DFDCSeed &seed, double theta_min, double theta_max)
504{
505 /// Find the z value of the vertex using the valid stereo hits from <i>seed</i>. The values
506 /// for phi_hit and pos.Z() are assumed to be valid as is the status of the
507 /// VALID_HIT bit in flags.
508 ///
509 /// This uses a histogramming technique that looks at the overlaps of the
510 /// z ranges subtended by each hit between the given theta limits.
511 /// The overlaps usually lead to a range of values for z_vertex. The limits
512 /// of these are stored in the z_min and z_max fields of the seed.
513 /// The centroid of the range is stored in the z_vertex field.
514
515 // We use a simple array to store our histogram here. We don't want to use
516 // ROOT histograms because they are not thread safe.
517 unsigned int Nbins = 300;
518 unsigned int hist[Nbins];
519 for(unsigned int i=0; i<Nbins; i++)hist[i] = 0; // clear histogram
520 double bin_width = 0.5; // bins are 5mm
521 double hist_low_limit = 0.0; // lower edge of histogram limits
522
523 // We effectively extend the theta_min and theta_max angles here
524 // a bit to include some error. The motivation is that if
525 // theta_min == theta_max that leads to z_min == z_max so there
526 // is little or no overlap of the z ranges of separate hits.
527 // For now, we hardwire this to 1 degree
528 double theta_err = 1.0/57.3;
529
530 // Loop over valid hits, filling the histogram
531 double r0 = seed.r0;
532 double tan_alpha_max = tan(theta_max + theta_err)/r0;
533 double tan_alpha_min = tan(theta_min - theta_err)/r0;
534 if(tan_alpha_min<0.0)tan_alpha_min=0.0;
535 for(unsigned int i=0; i<seed.hits.size(); i++){
536 DFDCTrkHit *trkhit = seed.hits[i];
537 if(!trkhit->flags&VALID_HIT)continue;
538 if(!trkhit->flags&ON_CIRCLE)continue;
539 if(!trkhit->flags&IN_THETA_RANGE)continue;
540
541 // Calculate upper and lower limits in z
542 double q_sign = seed.q>0.0 ? +1.0:-1.0;
543 double z_hit = trkhit->hit->pos.Z();
544 double zmin = z_hit - q_sign*trkhit->phi_hit/tan_alpha_min;
545 double zmax = z_hit - q_sign*trkhit->phi_hit/tan_alpha_max;
546 if(zmin>zmax){
547 double tmp = zmin;
548 zmin=zmax;
549 zmax=tmp;
550 }
551
552 // Copy z limits into trkhit
553 trkhit->zmin = zmin;
554 trkhit->zmax = zmax;
555
556 // Find index of bins corresponding to tmin and tmax
557 unsigned int imin = zmin<hist_low_limit ? 0:(unsigned int)floor((zmin-hist_low_limit)/bin_width);
558 unsigned int imax = zmax<hist_low_limit ? 0:(unsigned int)floor((zmax-hist_low_limit)/bin_width);
559
560 if(debug_level>3)_DBG_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
560<<" "<<" -- phi_hit="<<trkhit->phi_hit<<" z_hit="<<z_hit<<endl;
561 if(debug_level>3)_DBG_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
561<<" "<<" -- zmin="<<zmin<<" zmax="<<zmax<<endl;
562 if(debug_level>3)_DBG_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
562<<" "<<" -- imin="<<imin<<" imax="<<imax<<endl;
563
564 // If entire range of this hit is outside of the histogram limit
565 // then discard this hit.
566 if(imax<=0 || imin>=Nbins)continue;
567
568 // Clip limits of bin range to our histogram limits
569 if(imin<0)imin=0;
570 if(imin>=Nbins)imin=Nbins-1;
571 if(imax<0)imax=0;
572 if(imax>=Nbins)imax=Nbins-1;
573
574 // Increment all bins between imin and imax
575 for(unsigned int j=imin; j<=imax; j++)hist[j]++;
576 }
577
578 // Look for the indexes of the plateau
579 unsigned int istart=(unsigned int)((TARGET_Z_MIN-hist_low_limit)/bin_width);
580 unsigned int iend=0;
581 for(unsigned int i=1; i<Nbins; i++){
582
583 // Only look in Target area
584 double z = hist_low_limit + bin_width*(0.5+(double)i);
585 if(z<TARGET_Z_MIN || z>TARGET_Z_MAX)continue;
586
587 if(hist[i]>hist[istart]){
588 istart = i;
589 if(debug_level>3)_DBG_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
589<<" "<<" -- istart="<<istart<<" (z="<<hist_low_limit + bin_width*(0.5+(double)istart)<<" , N="<<hist[i]<<")"<<endl;
590 }
591 if(hist[i] == hist[istart])iend = i;
592 }
593
594 // If there are no entries in the histogram, then flag this seed as invalid
595 if(hist[istart]==0.0){
596 if(debug_level>3)_DBG_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
596<<" "<<" - No entries in z-vertex hist. Seed aborted."<<endl;
597 seed.valid=false;
598 }
599
600 // Calculate z limits
601 seed.z_min = hist_low_limit + bin_width*(0.5+(double)istart);
602 seed.z_max = hist_low_limit + bin_width*(0.5+(double)iend);
603 seed.z_vertex = (seed.z_max + seed.z_min)/2.0;
604 if(debug_level>3)_DBG_std::cerr<<"DTrackCandidate_factory_FDC.cc"<<":"<<
604<<" "<<"istart="<<istart<<" iend="<<iend<<" z_min="<<seed.z_min<<" z_max="<<seed.z_max<<" hits[istart]="<<hist[istart]<<endl;
605
606 // Mark all of the hits that have a theta range consistent with seed.theta
607 // and a z_vertex consistent with seed.z_vertex as being IN_Z_RANGE.
608 for(unsigned int i=0; i<seed.hits.size(); i++){
609 DFDCTrkHit *trkhit = seed.hits[i];
610 trkhit->flags &= ~IN_Z_RANGE;
611 if(!trkhit->flags&VALID_HIT)continue;
612 if(!trkhit->flags&ON_CIRCLE)continue;
613 if(!trkhit->flags&IN_THETA_RANGE)continue;
614 if(trkhit->zmin > seed.z_vertex)continue;
615 if(trkhit->zmax < seed.z_vertex)continue;
616 trkhit->flags |= IN_Z_RANGE;
617 }
618}
619