Bug Summary

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

Annotated Source Code

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