Bug Summary

File:libraries/TRACKING/DTrackCandidate_factory_CDC.cc
Location:line 1197, column 11
Description:Dereference of null pointer

Annotated Source Code

1// $Id$
2//
3// File: DTrackCandidate_factory_CDC.cc
4// Created: Thu Sep 6 14:47:48 EDT 2007
5// Creator: davidl (on Darwin Amelia.local 8.10.1 i386)
6//
7
8#include <algorithm>
9using namespace std;
10
11#include <DVector2.h>
12#include <CDC/DCDCTrackHit.h>
13#include "DHelicalFit.h"
14#include <HDGEOMETRY/DGeometry.h>
15#include <HDGEOMETRY/DMagneticFieldMap.h>
16
17#define BeamRMS1.0 1.0
18#define EPS1e-3 1e-3
19
20#include "DTrackCandidate_factory_CDC.h"
21
22#define TWO(c)(0x1u << (c)) (0x1u << (c))
23#define MASK(c)((unsigned int)(-1)) / ((0x1u << ((0x1u << (c))))
+ 1u) ((unsigned int)(-1)) / (TWO(TWO(c))(0x1u << ((0x1u << (c)))) + 1u)
24#define COUNT(x,c)((x) & ((unsigned int)(-1)) / ((0x1u << ((0x1u <<
(c)))) + 1u)) + (((x) >> ((0x1u << (c)))) & (
(unsigned int)(-1)) / ((0x1u << ((0x1u << (c)))) +
1u)) ((x) & MASK(c)((unsigned int)(-1)) / ((0x1u << ((0x1u << (c))))
+ 1u)) + (((x) >> (TWO(c)(0x1u << (c)))) & MASK(c)((unsigned int)(-1)) / ((0x1u << ((0x1u << (c))))
+ 1u))
25
26inline int bitcount (unsigned int n) {
27 n = COUNT(n, 0)((n) & ((unsigned int)(-1)) / ((0x1u << ((0x1u <<
(0)))) + 1u)) + (((n) >> ((0x1u << (0)))) & (
(unsigned int)(-1)) / ((0x1u << ((0x1u << (0)))) +
1u)) ;
28 n = COUNT(n, 1)((n) & ((unsigned int)(-1)) / ((0x1u << ((0x1u <<
(1)))) + 1u)) + (((n) >> ((0x1u << (1)))) & (
(unsigned int)(-1)) / ((0x1u << ((0x1u << (1)))) +
1u)) ;
29 n = COUNT(n, 2)((n) & ((unsigned int)(-1)) / ((0x1u << ((0x1u <<
(2)))) + 1u)) + (((n) >> ((0x1u << (2)))) & (
(unsigned int)(-1)) / ((0x1u << ((0x1u << (2)))) +
1u)) ;
30 n = COUNT(n, 3)((n) & ((unsigned int)(-1)) / ((0x1u << ((0x1u <<
(3)))) + 1u)) + (((n) >> ((0x1u << (3)))) & (
(unsigned int)(-1)) / ((0x1u << ((0x1u << (3)))) +
1u)) ;
31 n = COUNT(n, 4)((n) & ((unsigned int)(-1)) / ((0x1u << ((0x1u <<
(4)))) + 1u)) + (((n) >> ((0x1u << (4)))) & (
(unsigned int)(-1)) / ((0x1u << ((0x1u << (4)))) +
1u)) ;
32 /*n = COUNT(n, 5) ; for 64-bit integers */
33 return n ;
34}
35
36
37
38class DVector3_with_perp:public DVector3
39{
40 public:
41 DVector3_with_perp():DVector3(){CalcPerp();}
42 DVector3_with_perp(double x, double y, double z):DVector3(x,y,z){CalcPerp();}
43 double CalcPerp(void){
44 perp = Perp();
45 return perp;
46 }
47 double perp;
48};
49
50inline bool SortIntersections(const DVector3_with_perp &a,const DVector3_with_perp &b){
51 if (a.perp<b.perp) return true;
52 return false;
53}
54
55inline bool CDCSortByRdecreasing(DTrackCandidate_factory_CDC::DCDCTrkHit* const &hit1, DTrackCandidate_factory_CDC::DCDCTrkHit* const &hit2) {
56 // use the ring number to sort by R(decreasing) and then straw(increasing)
57 if(hit1->hit->wire->ring == hit2->hit->wire->ring){
58 return hit1->hit->wire->straw < hit2->hit->wire->straw;
59 }
60 return hit1->hit->wire->ring > hit2->hit->wire->ring;
61}
62inline bool CDCSortByRincreasing(DTrackCandidate_factory_CDC::DCDCTrkHit* const &hit1, DTrackCandidate_factory_CDC::DCDCTrkHit* const &hit2) {
63 // use the ring number to sort by R
64 return hit1->hit->wire->ring < hit2->hit->wire->ring;
65}
66inline bool CDCSortByZincreasing(DTrackCandidate_factory_CDC::DCDCTrkHit* const &hit1, DTrackCandidate_factory_CDC::DCDCTrkHit* const &hit2) {
67 // use the z_stereo position to sort
68 return hit1->z_stereo < hit2->z_stereo;
69}
70inline bool CDCSortByStereoPhiincreasing(DTrackCandidate_factory_CDC::DCDCTrkHit* const &hit1, DTrackCandidate_factory_CDC::DCDCTrkHit* const &hit2) {
71 // use the z_stereo position to sort
72 return hit1->phi_stereo < hit2->phi_stereo;
73}
74
75//------------------
76// init
77//------------------
78jerror_t DTrackCandidate_factory_CDC::init(void)
79{
80 MAX_ALLOWED_CDC_HITS = 10000;
81 MAX_SUBSEED_STRAW_DIFF = 1;
82 MIN_SEED_HITS = 2;
83 MAX_SUBSEED_LINKED_HITS = 12;
84 MAX_RING_SUBSEED_HITS = 4;
85 MAX_HIT_DIST = 4.0; // cm
86 MAX_SEED_TIME_DIFF = 1000.0; // ns
87 MAX_CDC_MATCH_ANGLE = 20.0; // degrees
88
89 MAX_SEED_LINK_ANGLE = M_PI3.14159265358979323846/6.0*57.3; // degrees
90 TARGET_Z_MIN = 50.0;
91 TARGET_Z_MAX = 80.0;
92 TARGET_Z=65.0;
93 DEBUG_LEVEL = 0;
94
95 FILTER_SEEDS=false;
96
97 // Initialize cdchits_by_superlayer with empty vectors for each superlayer
98 vector<DCDCTrkHit*> mt;
99 for(int i=0; i<5; i++)cdchits_by_superlayer.push_back(mt);
100
101 // Set the layer numbers defining the superlayer boundaries
102 superlayer_boundaries.push_back( 4);
103 superlayer_boundaries.push_back(12);
104 superlayer_boundaries.push_back(16);
105 superlayer_boundaries.push_back(24);
106 superlayer_boundaries.push_back(28);
107
108 return NOERROR;
109}
110
111//------------------
112// brun
113//------------------
114jerror_t DTrackCandidate_factory_CDC::brun(JEventLoop *eventLoop, int runnumber)
115{
116 DApplication* dapp=dynamic_cast<DApplication*>(eventLoop->GetJApplication());
117 bfield = dapp->GetBfield();
118
119 // get the geometry
120 const DGeometry *geom = dapp->GetDGeometry(runnumber);
121
122 geom->GetTargetZ(TARGET_Z);
123 double zrange;
124 geom->GetTargetLength(zrange);
125 TARGET_Z_MIN=TARGET_Z-0.5*zrange;
126 TARGET_Z_MAX=TARGET_Z+0.5*zrange;
127
128
129 gPARMS->SetDefaultParameter("TRKFIND:MAX_ALLOWED_CDC_HITS", MAX_ALLOWED_CDC_HITS);
130 gPARMS->SetDefaultParameter("TRKFIND:MAX_SUBSEED_STRAW_DIFF", MAX_SUBSEED_STRAW_DIFF);
131 gPARMS->SetDefaultParameter("TRKFIND:MIN_SEED_HITS", MIN_SEED_HITS);
132 gPARMS->SetDefaultParameter("TRKFIND:MAX_SUBSEED_LINKED_HITS", MAX_SUBSEED_LINKED_HITS);
133 gPARMS->SetDefaultParameter("TRKFIND:MAX_RING_SUBSEED_HITS", MAX_RING_SUBSEED_HITS);
134 gPARMS->SetDefaultParameter("TRKFIND:MAX_HIT_DIST", MAX_HIT_DIST);
135 gPARMS->SetDefaultParameter("TRKFIND:MAX_SEED_TIME_DIFF", MAX_SEED_TIME_DIFF);
136 gPARMS->SetDefaultParameter("TRKFIND:MAX_CDC_MATCH_ANGLE", MAX_CDC_MATCH_ANGLE);
137
138 gPARMS->SetDefaultParameter("TRKFIND:MAX_SEED_LINK_ANGLE", MAX_SEED_LINK_ANGLE);
139 gPARMS->SetDefaultParameter("TRKFIND:DEBUG_LEVEL", DEBUG_LEVEL);
140 gPARMS->SetDefaultParameter("TRKFIND:FILTER_SEEDS",FILTER_SEEDS);
141
142 MAX_HIT_DIST2 = MAX_HIT_DIST*MAX_HIT_DIST;
143
144 return NOERROR;
145}
146
147//------------------
148// erun
149//------------------
150jerror_t DTrackCandidate_factory_CDC::erun(void)
151{
152 return NOERROR;
153}
154
155//------------------
156// fini
157//------------------
158jerror_t DTrackCandidate_factory_CDC::fini(void)
159{
160 for(unsigned int i=0; i<cdctrkhits.size(); i++)delete cdctrkhits[i];
161 return NOERROR;
162}
163
164//------------------
165// evnt
166//------------------
167jerror_t DTrackCandidate_factory_CDC::evnt(JEventLoop *loop, int eventnumber)
168{
169 // Get CDC hits
170 unsigned int numcdchits=0;
171 if (GetCDCHits(loop,numcdchits)!=NOERROR){
172 return RESOURCE_UNAVAILABLE;
173 }
174
175 // Find all seeds in all 3 axial superlayers
176 vector<DCDCSeed> seeds_sl5;
177 vector<DCDCSeed> seeds_sl3;
178 vector<DCDCSeed> seeds_sl1;
179 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
179<<" "<<"========= SL5 =========="<<endl;
180 FindSeeds(cdchits_by_superlayer[5-1], seeds_sl5);
181 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
181<<" "<<"========= SL3 =========="<<endl;
182 FindSeeds(cdchits_by_superlayer[3-1], seeds_sl3);
183 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
183<<" "<<"========= SL1 =========="<<endl;
184 FindSeeds(cdchits_by_superlayer[1-1], seeds_sl1);
185
186 // Link seeds using average phi of seed hits
187 vector<DCDCSeed> seeds_tmp, seeds_tmpSL5,seeds_tmpSL3_SL5,seeds;
188 LinkSeeds(seeds_sl5, seeds_sl3, seeds_tmp, MAX_SUBSEED_LINKED_HITS);
189
190 // Put aside the seeds that only have hits in SL5
191 for (unsigned int i=0;i<seeds_tmp.size();i++){
192 unsigned int max_ind=seeds_tmp[i].hits.size()-1;
193 if (seeds_tmp[i].hits[max_ind]->hit->wire->ring>24){
194 seeds_tmpSL5.push_back(seeds_tmp[i]);
195 }
196 else{
197 seeds_tmpSL3_SL5.push_back(seeds_tmp[i]);
198 }
199 }
200
201 // Link seeds between SL1 and SL3+SL5
202 LinkSeeds(seeds_tmpSL3_SL5, seeds_sl1, seeds, 2*MAX_SUBSEED_LINKED_HITS);
203
204 // Insert the SL5 seeds into the full list of seeds
205 seeds.insert(seeds.begin(),seeds_tmpSL5.begin(),seeds_tmpSL5.end());
206
207 // Check to add lone hits in SL3 to seeds with only SL1 hits
208 PickupUnmatched(seeds);
209
210 // Drop seeds containing hits from only a single axial superlayer besides SL1
211 DropIncompleteSeeds(seeds);
212
213 if (seeds.size()>0){
214 // Create vectors to store bit pattern of hits used by the seeds
215 for (unsigned int i=0;i<seeds.size();i++){
216 seeds[i].HitBitPattern.clear();
217 seeds[i].HitBitPattern.resize(numcdchits/(8*sizeof(unsigned int))+1);
218 }
219
220 // Fit linked seeds to circles
221 for(unsigned int i=0; i<seeds.size(); i++){
222 if(DEBUG_LEVEL>5)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
222<<" "<<"-- Starting fit for seed "<<i<<" Nhits="<<seeds[i].hits.size()<<" phi_avg="<<seeds[i].phi_avg<<endl;
223 seeds[i].valid = FitCircle(seeds[i]);
224 }
225
226 // Filter out duplicates of seeds by clearing their "valid" flags
227 FilterCloneSeeds(seeds);
228
229 // Extend seeds into stereo layers and do fit to find z and theta
230 for(unsigned int i=0; i<seeds.size(); i++){
231 DCDCSeed &seed = seeds[i];
232 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
232<<" "<<"----- Seed "<<i<<" ------"<<endl;
233 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
233<<" "<<"seed.fit.phi="<<seed.fit.phi<<endl;
234 if(!seed.valid)continue;
235
236 // Add stereo hits to seed
237 AddStereoHits(cdchits_by_superlayer[2-1], seed);
238
239 // If there are hits in the first axial layers but no hits in the first
240 // stereo layers, then this is unlikely to be a valid (fittable)
241 // candidate...
242 unsigned int first_stereo_count=seed.stereo_hits.size();
243 if (first_stereo_count==0 && seed.hits[0]->hit->wire->ring<5){
244 seed.valid=false;
245 continue;
246 }
247
248 // Go on to the next set of stereo layers
249 AddStereoHits(cdchits_by_superlayer[4-1], seed);
250
251 // If no stereo hits were found for this seed, then
252 // we can't fit it.
253 if(seed.stereo_hits.size()==0){
254 seed.valid=false;
255 continue;
256 }
257
258 if (FindThetaZRegression(seed)!=NOERROR){
259 // If the linear regression doesn't work try the histogramming method
260 // Fit stereo hits to get theta and vertex z position
261 FindThetaZ(seed);
262 if(!seed.valid){
263 //continue;
264
265 // Assume that the track came from one end or the other
266 // of the target and use a point in one of the stereo
267 // layers to estimate tanl
268 if (seed.z_vertex>TARGET_Z_MAX)
269 seed.z_vertex=TARGET_Z_MAX;
270 else
271 seed.z_vertex=TARGET_Z_MIN;
272 double x=seed.stereo_hits[0].x_stereo;
273 double y=seed.stereo_hits[0].y_stereo;
274 double ratio=sqrt(x*x+y*y)/2./seed.fit.r0;
275 if (ratio<1.){
276 double tanl=(seed.stereo_hits[0].z_stereo-seed.z_vertex)/
277 (2.*seed.fit.r0*asin(ratio));
278 seed.theta=M_PI_21.57079632679489661923-atan(tanl);
279 }
280
281 }
282 }
283 }
284
285 if (FILTER_SEEDS){
286 //Look for seeds that share hits and mark a seed as invalid if it shares more
287 //than a certain fraction (to be determined) of hits with another seed and the
288 //chisq/df from the circle fit is worse than that of the other fit by a
289 //certain amount
290 unsigned int numwords=seeds[0].HitBitPattern.size();
291 for(unsigned int i=0; i<seeds.size()-1; i++){
292 if (seeds[i].valid){
293 for (unsigned int j=i+1;j<seeds.size();j++){
294 if (seeds[j].valid){
295 int numcommon=0;
296 for (unsigned int k=0;k<numwords;k++){
297 numcommon
298 +=bitcount(seeds[i].HitBitPattern[k]&seeds[j].HitBitPattern[k]);
299 }
300 double hitratio=double(numcommon)/
301 double(seeds[i].hits.size()+seeds[i].stereo_hits.size());
302 double chisq_nu_1=seeds[i].fit.chisq/seeds[i].fit.ndof;
303 double chisq_nu_2=seeds[j].fit.chisq/seeds[j].fit.ndof;
304 if (hitratio>0.9 || (hitratio>0.5 && fabs(chisq_nu_1-chisq_nu_2)>1)){
305 if (chisq_nu_1 > chisq_nu_2) seeds[i].valid=false;
306 else seeds[j].valid=false;
307 }
308 }
309 }
310 }
311 }
312 }
313
314 // Put the good seeds in the list of cdc track candidates
315 for(unsigned int i=0; i<seeds.size(); i++){
316 DCDCSeed seed = seeds[i];
317 if (seed.valid){
318 // The following is from a fit of ratio of thrown to reconstructed
319 // transverse momentum vs. theta for the 1400A field
320 //double par[] = {0.984463, 0.150759, -0.414933, 0.257472, -0.055801};
321 //double theta = seed.theta;
322 //double ff = par[0]+theta*(par[1]+theta*(par[2]+theta*(par[3]+theta*par[4])));
323
324 DVector3 pos, mom;
325 GetPositionAndMomentum(seed,pos,mom);
326
327 //Make a track candidate from results
328 DTrackCandidate *can = new DTrackCandidate;
329
330 can->setPosition(pos);
331 can->setMomentum(mom);
332 //can->setCharge(seed.q);
333 can->setCharge(seed.fit.q);
334
335 for (unsigned int n=0;n<seed.hits.size();n++){
336 const DCDCTrackHit *cdchit=(seed.hits[n])->hit;
337 can->used_cdc_indexes.push_back(seed.hits[n]->index);
338 can->AddAssociatedObject(cdchit);
339 }
340 for (unsigned int n=0;n<seed.stereo_hits.size();n++){
341 const DCDCTrackHit *cdchit=(seed.stereo_hits[n]).hit;
342 can->used_cdc_indexes.push_back((seed.stereo_hits[n]).index);
343 can->AddAssociatedObject(cdchit);
344 }
345
346
347 //can->q = can->charge();
348 //can->phi = mom.Phi();
349 //can->theta = mom.Theta();
350 //can->p_trans = mom.Pt();
351 //can->p = mom.Mag();
352 //can->z_vertex = pos.Z();
353 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
353<<" "<<"Final Candidate parameters: p="<<mom.Mag()<<" theta="<<mom.Theta()<<" phi="<<mom.Phi()<<" z="<<pos.Z()<<endl;
354
355 _data.push_back(can);
356 }
357 }
358 }
359
360 return NOERROR;
361}
362
363//------------------
364// GetCDCHits
365//------------------
366jerror_t DTrackCandidate_factory_CDC::GetCDCHits(JEventLoop *loop,
367 unsigned int &numhits)
368{
369 // Delete old hits
370 for(unsigned int i=0; i<cdctrkhits.size(); i++)delete cdctrkhits[i];
371 cdctrkhits.clear();
372 for(unsigned int i=0; i<cdchits_by_superlayer.size(); i++)cdchits_by_superlayer[i].clear();
373
374
375 // Get the "raw" hits. These already have the wire associated with them.
376 vector<const DCDCTrackHit*> cdctrackhits;
377 loop->Get(cdctrackhits);
378 numhits=cdctrackhits.size();
379
380 // If there are no hits, then bail now
381 if(cdctrackhits.size()==0)return RESOURCE_UNAVAILABLE;
382
383 // If there are too many hits, bail with a warning message
384 if(cdctrackhits.size()>MAX_ALLOWED_CDC_HITS){
385 _DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
385<<" "<<"Too many hits in CDC ("<<cdctrackhits.size()<<", max="<<MAX_ALLOWED_CDC_HITS<<")! Track finding in CDC bypassed for event "<<loop->GetJEvent().GetEventNumber()<<endl;
386 cdctrackhits.clear();
387 return UNRECOVERABLE_ERROR;
388 }
389
390 // Create DCDCTrkHit objects out of these.
391 int oldwire = -1;
392 for(unsigned int i=0; i<cdctrackhits.size(); i++){
393 // Add to "master" list
394 // ONLY FIRST HIT OF A WIRE
395 int newwire = cdctrackhits[i]->wire->ring*1000 + cdctrackhits[i]->wire->straw;
396 if (newwire!=oldwire){
397
398 oldwire = newwire;
399
400 // Add to "master" list
401 DCDCTrkHit *cdctrkhit = new DCDCTrkHit;
402 cdctrkhit->index=i;
403 cdctrkhit->hit = cdctrackhits[i];
404 cdctrkhit->flags = cdctrkhit->hit->is_stereo==false ? NONE:IS_STEREO;
405 cdctrkhit->flags |= NOISE; // (see below)
406
407 cdctrkhits.push_back(cdctrkhit);
408
409 // Sort into list of hits by superlayer
410#if 1
411 for(unsigned int j=0; j<superlayer_boundaries.size(); j++){
412 if(cdctrkhit->hit->wire->ring<=superlayer_boundaries[j]){
413 cdchits_by_superlayer[j].push_back(cdctrkhit);
414 break;
415 }
416 }
417#else
418 if(cdctrkhit->hit->wire->ring<=4)cdchits_by_superlayer[0].push_back(cdctrkhit);
419 else if(cdctrkhit->hit->wire->ring<= 12)cdchits_by_superlayer[1].push_back(cdctrkhit);
420 else if(cdctrkhit->hit->wire->ring<=16)cdchits_by_superlayer[2].push_back(cdctrkhit);
421 else if(cdctrkhit->hit->wire->ring<=24)cdchits_by_superlayer[3].push_back(cdctrkhit);
422 else if(cdctrkhit->hit->wire->ring<=28)cdchits_by_superlayer[4].push_back(cdctrkhit);
423#endif
424 }
425 }
426
427 // Sort the individual superlayer lists by decreasing values of R
428 for(unsigned int i=0; i<cdchits_by_superlayer.size(); i++){
429 sort(cdchits_by_superlayer[i].begin(), cdchits_by_superlayer[i].end(), CDCSortByRdecreasing);
430 }
431
432 // Filter out noise hits. All hits are initially flagged as "noise".
433 // Hits with a neighbor within MAX_HIT_DIST have their noise flags cleared.
434 for(unsigned int i=0; i<cdctrkhits.size(); i++){ // cut above should ensure cdctrkhits.size() is at least 1
435 DCDCTrkHit *trkhit1 = cdctrkhits[i];
436 if(!(trkhit1->flags & NOISE))continue; // this hit already not marked for noise
437 for(unsigned int j=0; j<cdctrkhits.size(); j++){
438 if(j==i)continue;
439 double d2 = trkhit1->Dist2(cdctrkhits[j]);
440 if(d2<9.0*MAX_HIT_DIST2){
441 trkhit1->flags &= ~NOISE;
442 cdctrkhits[j]->flags &= ~NOISE;
443 break;
444 }// if
445 }// j
446 }// i
447 return NOERROR;
448}
449
450//------------------
451// FindSeeds
452//------------------
453void DTrackCandidate_factory_CDC::FindSeeds(vector<DCDCTrkHit*> &hits, vector<DCDCSeed> &seeds)
454{
455 // Sort through hits ring by ring to find sub-seeds from neighboring wires
456 // in the same ring. What we want is a list of sub-seeds for each ring.
457 // Each sub-seed is a list of adjacent (or nearly adjacent) hits. Note
458 // that hits are ordered by ring, then straw. We ignore the boundary in
459 // straw number for now which will cause some sub-seeds to be listed as
460 // 2 separate ones. These will be recombined later.
461 DCDCSeed subseed;
462 vector<DCDCSeed > subseeds;
463 vector<vector<DCDCSeed > > rings;
464 int last_ring = -1;
465 for(unsigned int i=0; i<hits.size(); i++){
466 DCDCTrkHit *trkhit = hits[i];
467
468 // Clear USED, IN_SEED, IN_LINE, IN_TRACK flags for all these hits
469 trkhit->flags &= ~(USED|IN_SEED|IN_LINE|IN_TRACK);
470
471 //if(trkhit->flags & NOISE)continue;
472
473 // Check if ring number has changed.
474 if(trkhit->hit->wire->ring != last_ring){
475 if(subseed.hits.size()!=0)subseeds.push_back(subseed);
476 if(subseeds.size()!=0)rings.push_back(subseeds);
477 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
477<<" "<<" subseed hits:"<<subseed.hits.size()<<" subseeds:"<<subseeds.size()<<endl;
478 subseeds.clear();
479 subseed.hits.clear();
480 subseed.hits.push_back(trkhit);
481 subseed.linked=false;
482 last_ring = trkhit->hit->wire->ring;
483 continue;
484 }
485
486 // Check if this hit is a neighbor of the last hit added to the subseed
487 if((unsigned int)abs(subseed.hits[subseed.hits.size()-1]->hit->wire->straw - trkhit->hit->wire->straw)>MAX_SUBSEED_STRAW_DIFF){
488 if(subseed.hits.size()!=0)subseeds.push_back(subseed);
489 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
489<<" "<<" subseed hits:"<<subseed.hits.size()<<endl;
490 subseed.hits.clear();
491 subseed.linked=false;
492 }
493 subseed.hits.push_back(trkhit);
494 }
495 if(subseed.hits.size()!=0)subseeds.push_back(subseed);
496 if(subseeds.size()!=0)rings.push_back(subseeds);
497 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
497<<" "<<" subseed hits:"<<subseed.hits.size()<<" subseeds:"<<subseeds.size()<<endl;
498 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
498<<" "<<"rings: "<<rings.size()<<endl;
499
500 // If we have no "rings", then there must be no seeds. Bail now
501 if(rings.size()==0)return;
502
503 // Loop over rings
504 for(ringiter ring =rings.begin(); ring!=rings.end(); ring++){
505 vector<DCDCSeed> &subseeds = *ring;
506 ringiter next_ring = ring;
507 next_ring++;
508
509 // Loop over subseeds of this ring
510 for(unsigned int j=0; j<subseeds.size(); j++){
511 if(subseeds[j].linked)continue;
512 if(subseeds[j].hits.size()>MAX_RING_SUBSEED_HITS){
513 if(DEBUG_LEVEL>4)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
513<<" "<<"rejecting subseed for having too many hits in a single layer ("<<subseeds[j].hits.size()<<" max="<<MAX_RING_SUBSEED_HITS<<")"<<endl;
514 continue;
515 }
516
517 // This subseed was never used. Start a new seed
518 vector<DCDCSeed*> parent;
519 parent.push_back(&subseeds[j]);
520 LinkSubSeeds(parent, next_ring, rings.end(), seeds);
521 }
522 }
523
524 // Mark all seeds as valid
525 for(unsigned int i=0; i<seeds.size(); i++){
526 seeds[i].valid = true;
527 }
528}
529
530//------------------
531// LinkSubSeeds
532//------------------
533void DTrackCandidate_factory_CDC::LinkSubSeeds(vector<DCDCSeed*> &parent, ringiter ring, ringiter ringend, vector<DCDCSeed> &seeds)
534{
535 /// Combine subseeds from layers in a single superlayer into seeds.
536 ///
537 /// This a a re-entrant routine (i.e. it calls itself recursively). Upon
538 /// entry, <i>parent</i> contains a list of pointers to all of the subseeds
539 /// from the rings outside of <i>ring</i> that are to be combined into
540 /// a seed. This will search through all subseeds of <i>ring</i> and if
541 /// any are found that can extend the parent, a copy of parent is made,
542 /// the current subseed of this ring is added to it, and then it is
543 /// passed on to another call to this routine. If no matches are found
544 /// (which will necessarily be the case for the inner most layer), then
545 /// the subseeds in <i>parent</i> will be combined into a single seed and
546 /// that added to <i>seeds</i>.
547
548 // Make sure parent has at least one subseed
549 if(parent.size()==0){
550 _DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
550<<" "<<"parent has no subseeds!!"<<endl;
551 return;
552 }
553
554 // Set flag to keep track of whether this is the end of the seed or not
555 bool seed_extended = false;
556 if(ring!=ringend){
557
558 // Last subseed in parent list is the one we need to compare to
559 DCDCSeed *parent_subseed = parent[parent.size()-1];
560 double r_parent = parent_subseed->hits[0]->hit->wire->origin.Perp();
561
562 // Loop over subseeds in this ring
563 vector<DCDCSeed> &subseeds = (*ring);
564 ring++; // increment ring iterator to point to next level down in case we recall ouself below
565
566 for(unsigned int i=0; i<subseeds.size(); i++){
567 // Find the difference between this subseed's distance from the beamline
568 // and the parent's.
569 double dr = r_parent - subseeds[i].hits[0]->hit->wire->origin.Perp();
570
571 // Check if this subseed is close enough to the parent's to extend it
572 if(parent_subseed->MinDist2(subseeds[i])<(dr*dr+MAX_HIT_DIST2)){
573 vector<DCDCSeed*> myparent = parent;
574 myparent.push_back(&subseeds[i]);
575 subseeds[i].linked = true;
576 LinkSubSeeds(myparent, ring, ringend, seeds);
577 seed_extended = true;
578 }
579 }
580 }
581
582 // Check if this is the end of the line.
583 if(!seed_extended){
584 // This is the end of this seed. Combine the subseeds into a single
585 // seed and add it to the list.
586 DCDCSeed seed;
587 seed.hits.clear();
588 for(unsigned int i=0; i<parent.size(); i++){
589 seed.Merge(*parent[i]);
590 }
591 if(seed.hits.size()>=MIN_SEED_HITS){
592 seeds.push_back(seed);
593 }else{
594 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
594<<" "<<"rejecting seed due to too few hits (have "<<seed.hits.size()<<" need "<<MIN_SEED_HITS<<")"<<endl;
595 }
596 }
597}
598
599//------------------
600// LinkSeeds
601//------------------
602void DTrackCandidate_factory_CDC::LinkSeeds(vector<DCDCSeed> &in_seeds1, vector<DCDCSeed> &in_seeds2, vector<DCDCSeed> &seeds, unsigned int max_linked_hits)
603{
604 /// Loop over the two input sets of seeds and compare the phi angles of
605 /// their first and last hits. The closest combination is checked to see
606 /// if we should combine them into a new a new seed. Any seeds from either
607 /// list that are not linked, but have enough hits to be a full seed in and
608 /// of themselves will be added to the <i>seeds</i> list.
609 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
609<<" "<<"Linking seeds: Num seeds in list 1:"<<in_seeds1.size()<<" Num seeds in list 2:"<<in_seeds2.size()<<endl;
610
611 // Clear all "linked" flags
612 for(unsigned int i=0; i< in_seeds1.size(); i++)in_seeds1[i].linked=false;
613 for(unsigned int i=0; i< in_seeds2.size(); i++)in_seeds2[i].linked=false;
614
615 for(unsigned int i=0; i< in_seeds1.size(); i++){
616 vector<DCDCTrkHit*> &hits1 = in_seeds1[i].hits;
617 if(DEBUG_LEVEL>5)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
617<<" "<<"hits1.size()="<<hits1.size()<<endl;
618 if(hits1.size()<1)continue;
619
620 for(unsigned int j=0; j< in_seeds2.size(); j++){
621 vector<DCDCTrkHit*> &hits2 = in_seeds2[j].hits;
622 if(DEBUG_LEVEL>5)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
622<<" "<<" hits2.size()="<<hits2.size()<<endl;
623 if(hits2.size()<1)continue;
624
625 // Find the hits from the seeds that are closest in R
626 const DCDCWire *wire1a = hits1[0]->hit->wire;
627 const DCDCWire *wire1b = hits1[hits1.size()-1]->hit->wire;
628 const DCDCWire *wire2a = hits2[0]->hit->wire;
629 const DCDCWire *wire2b = hits2[hits2.size()-1]->hit->wire;
630 const DVector3 *pos1, *pos2;
631 if(wire1a->ring > wire2a->ring){
632 // seed1 is outside seed2
633 pos1 = &(wire1a->ring > wire1b->ring ? wire1b:wire1a)->origin;
634 pos2 = &(wire2a->ring < wire2b->ring ? wire2b:wire2a)->origin;
635 }else{
636 // seed2 is outside seed1
637 pos1 = &(wire1a->ring < wire1b->ring ? wire1b:wire1a)->origin;
638 pos2 = &(wire2a->ring > wire2b->ring ? wire2b:wire2a)->origin;
639 }
640
641 // Compare the phi values of the first and last points of the two seeds
642 // to see if they are close enough to link.
643 double dphi = fabs(pos1->Phi() - pos2->Phi());
644 while(dphi>M_PI3.14159265358979323846)dphi-=2.0*M_PI3.14159265358979323846;
645 if(fabs(dphi*57.3)<MAX_SEED_LINK_ANGLE){
646
647 DCDCSeed seed = in_seeds1[i];
648 seed.Merge(in_seeds2[j]);
649 seeds.push_back(seed);
650
651 // Mark all of the hits in the new, merged seed as USED. Some will
652 // already have the flag set by having been previously merged. The
653 // USED flag is used later when matching SL1 only seeds with stray
654 // hits in SL3.
655 for(unsigned int k=0; k<seed.hits.size(); k++){
656 seed.hits[k]->flags |= USED;
657 }
658
659 // OK, at this point we have linked the seeds and we *may* want to
660 // to set thier "linked" flags. The linked flags are really used
661 // to decide later if the partial seed should be promoted to
662 // a full seed. The idea here is to only flag the seed as "linked"
663 // if the link is strong, otherwise leaved it marked as unlinked
664 // so it has the option of becoming another track candidate at
665 // the next level. Note that even if we set a linked flag here,
666 // the partial seed can still be linked with other partials in
667 // this loop.
668 //
669 // The criteria for a strong link are:
670 // 1. the number of hits in the partner's partial seed is not
671 // greater than max_linked_hits (which is passed to us as an
672 // argument)
673 // 2. The minimum distance between the partial seeds' ends is less
674 // than 3*MAX_HIT_DIST in the phi direction.
675
676 // Find distance between seeds' ends
677 double dr = fabs(pos1->Perp() - pos2->Perp());
678 DVector3 d = *pos1 - *pos2;
679 double dist_phi2 = d.Perp2() - dr*dr;
680
681 if(dist_phi2 < 9.0*MAX_HIT_DIST2){
682
683 // Mark the seeds has having been linked. Here we do something
684 // a little funny. We only mark a seed as linked if the other
685 // seed doesn't have too many hits. This is because we can get
686 // "super seeds" from a single superlayer made from many hits,
687 // often due to crossing tracks.
688 if(hits1.size()<=max_linked_hits)in_seeds2[j].linked = true;
689 if(hits2.size()<=max_linked_hits)in_seeds1[i].linked = true;
690 }else{
691 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
691<<" "<<" not marking seeds as linked (dist_phi="<<sqrt(dist_phi2)<<" cm > "<<sqrt(9.0*MAX_HIT_DIST2)<<" cm)"<<endl;
692 }
693 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
693<<" "<<" linked seeds "<<i<<" ("<<hits1.size()<<" hits) and "<<j<<" ("<<hits2.size()<<" hits) dist_phi="<<sqrt(dist_phi2)<<endl;
694 }else{
695 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
695<<" "<<" rejected link between "<<i<<" and "<<j<<" due to avg. phi (fabs("<<in_seeds1[i].phi_avg*57.3<<" - "<<in_seeds2[j].phi_avg*57.3<<")>="<<MAX_SEED_LINK_ANGLE<<")"<<endl;
696 }
697 }
698 }
699
700 // Any seeds that weren't linked from either list should be promoted
701 // to full seeds.
702 for(unsigned int i=0; i< in_seeds1.size(); i++){
703 if(!in_seeds1[i].linked){
704 seeds.push_back(in_seeds1[i]);
705 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
705<<" "<<"Promoting seed "<<i<<" from list 1"<<endl;
706 }
707 }
708 for(unsigned int i=0; i< in_seeds2.size(); i++){
709 if(!in_seeds2[i].linked){
710 seeds.push_back(in_seeds2[i]);
711 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
711<<" "<<"Promoting seed "<<i<<" from list 2"<<endl;
712 }
713 }
714
715 // Calculate average time for all axial hits in each seed and set
716 // all seeds to a charge of +1
717 for(unsigned int i=0; i< seeds.size(); i++){
718 DCDCSeed &seed = seeds[i];
719 seed.tdrift_avg = 0.0;
720 unsigned int Nin_time=0;
721 for(unsigned int j=0; j<seed.hits.size(); j++){
722 double tdrift = seed.hits[j]->hit->tdrift;
723 if(tdrift>1000.0){
724 seed.hits[j]->flags |= OUT_OF_TIME;
725 continue;
726 }
727 seed.tdrift_avg += tdrift;
728 Nin_time++;
729 }
730 seed.tdrift_avg /= (double)Nin_time;
731 }
732}
733
734//------------------
735// FitCircle
736//------------------
737bool DTrackCandidate_factory_CDC::FitCircle(DCDCSeed &seed)
738{
739 /// Do a quick fit of the 2-D projection of the axial hits
740 /// in seed (seed should contain *only* axial hits) to a
741 /// circle. Determine the sign of the charge (and correspondingly
742 /// the initial phi angle) assuming that
743 /// the majority of hits come from the outgoing part of the
744 /// track. If the resulting circle passes within
745 /// MAX_HIT_DIST of more the majority of hits, then true
746 /// is returned indicating success. Otherwise, false is
747 /// returned indicating failure and that the seed should be
748 /// discarded.
749
750 // Loop over hits in seed and add them to the seed's DHelicalFit object
751 seed.fit.Clear();
752 for(unsigned int j=0; j<seed.hits.size(); j++){
753 unsigned int numbits=8*sizeof(unsigned int);
754 seed.HitBitPattern[seed.hits[j]->index/numbits]
755 |=0x1u<<seed.hits[j]->index%numbits;
756
757 if(seed.hits[j]->flags&OUT_OF_TIME)continue;
758
759 const DVector3 &pos = seed.hits[j]->hit->wire->origin;
760 seed.fit.AddHitXYZ(pos.X(), pos.Y(),pos.Z());
761 }
762
763 // Try and fit the circle using a Riemann fit. If
764 // it fails, try a basic fit with QuickFit.
765 if(seed.fit.FitCircleRiemann(TARGET_Z,BeamRMS1.0)!=NOERROR
766 /* || seed.fit.chisq>20*/
767 ){
768 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
768<<" "<<"Riemann fit failed. Attempting regression fit..."<<endl;
769 if(seed.fit.FitCircle()!=NOERROR
770 /*|| seed.fit.chisq>20*/
771 ){
772 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
772<<" "<<"Regression circle fit failed. Trying straight line."<<endl;
773 if(seed.fit.FitCircleStraightTrack()!=NOERROR){
774 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
774<<" "<<"Trying straight line fit failed. Giving up."<<endl;
775 return false;
776 }
777 }else{
778 seed.fit.GuessChargeFromCircleFit(); // for Riemann fit
779 }
780 }else{
781 seed.fit.GuessChargeFromCircleFit(); // for regression fit
782 }
783
784 // Check if majority of hits are close to circle.
785 double x0 = seed.fit.x0;
786 double y0 = seed.fit.y0;
787 double r0 = seed.fit.r0;
788
789 unsigned int N=0;
790 for(unsigned int i=0; i<seed.hits.size(); i++){
791 if(seed.hits[i]->flags&OUT_OF_TIME){
792 if(DEBUG_LEVEL>6)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
792<<" "<<"discarding out of time hit"<<endl;
793 continue;
794 }
795 double dx = seed.hits[i]->hit->wire->origin.X() - x0;
796 double dy = seed.hits[i]->hit->wire->origin.Y() - y0;
797 double d = sqrt(dx*dx + dy*dy);
798 if(DEBUG_LEVEL>15)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
798<<" "<<"dist="<<d-r0<<endl;
799 if(fabs(d-r0)<=MAX_HIT_DIST)N++;
800 }
801
802 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
802<<" "<<"Circle fit: Nhits="<<seed.fit.GetNhits()<<" p_trans="<<seed.fit.p_trans<<" q="<<seed.fit.q<<" N="<<N<<" phi="<<seed.fit.phi<<endl;
803 if(N<MIN_SEED_HITS){
804 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
804<<" "<<"Rejected circle fit due to too few hits on track (N="<<N<<" MIN_SEED_HITS="<<MIN_SEED_HITS<<")"<<endl;
805 return false;
806 }
807
808 if(N<seed.hits.size()/2){
809 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
809<<" "<<"Rejected circle fit due to minority of hits on track (N="<<N<<" seed.hits.size()/2="<<seed.hits.size()/2<<")"<<endl;
810 return false;
811 }
812 return true;
813}
814
815//------------------
816// PickupUnmatched
817//------------------
818void DTrackCandidate_factory_CDC::PickupUnmatched(vector<DCDCSeed> &seeds)
819{
820 /// Look for single hits in superlayers that did not have enough
821 /// neighbors to make seeds of their own, but could be part of
822 /// a seed made in another axial superlayer. This is mainly here
823 /// to address tracks that have a single hit in SL3 that is
824 /// otherwise discarded. Seeds with only hits in SL1 often
825 /// don't enough information to form a good candidate so we try
826 /// and pick up these single hits in order to get a good handle
827 /// on the candidate.
828
829 // Loop through seeds, looking for ones with only SL1 hits
830 for(unsigned int i=0; i<seeds.size(); i++){
831 if(!seeds[i].valid)continue;
832 vector<DCDCTrkHit*> &hits = seeds[i].hits;
833 if(hits.size()<1)continue;
834 bool has_non_SL1_hit = false;
835 for(unsigned int j=0; j<hits.size(); j++){
836 if(hits[j]->hit->wire->ring>superlayer_boundaries[0]){
837 has_non_SL1_hit = true;
838 break;
839 }
840 }
841 if(has_non_SL1_hit)continue;
842 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
842<<" "<<"seed "<<i<<" has only SL1 hits. Looking for stray hits to match it with ..."<<endl;
843
844 // This seed must have only SL1 hits, look for unused SL3 hits
845 // at a phi angle consistent with the outermost hit.
846 bool added_hit = false;
847 const DCDCWire *sl1_wire = hits[hits.size()-1]->hit->wire;
848 for(unsigned int j=0; j<cdctrkhits.size(); j++){
849 DCDCTrkHit *sl3_hit = cdctrkhits[j];
850 if(sl3_hit->flags&USED)continue;
851 if(sl3_hit->hit->wire->ring<=superlayer_boundaries[1])continue;
852 if(sl3_hit->hit->wire->ring>superlayer_boundaries[2])continue;
853 double a = sl1_wire->origin.Angle(sl3_hit->hit->wire->origin);
854 if(fabs(a)<MAX_CDC_MATCH_ANGLE/57.3){
855 sl3_hit->flags |= USED;
856 hits.push_back(sl3_hit);
857 added_hit = true;
858 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
858<<" "<<"Adding stray SL3 hit to SL1 seed a="<<a*57.3<<" flags="<<sl3_hit->flags<<" ring="<<sl3_hit->hit->wire->ring<<endl;
859 }
860 }
861
862 if(added_hit){
863 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
863<<" "<<" found SL3 CDC hit to add to seed "<<i<<endl;
864 continue; // We found at least one SL3 hit that was added to the seed
865 }
866
867 }
868}
869
870//------------------
871// DropIncompleteSeeds
872//------------------
873void DTrackCandidate_factory_CDC::DropIncompleteSeeds(vector<DCDCSeed> &seeds)
874{
875 /// Look for seeds that have hits only in SL3 or only in SL5 and drop them.
876
877 int iSL1_lo = 1;
878 int iSL1_hi = superlayer_boundaries[0];
879 int iSL3_lo = superlayer_boundaries[1]+1;
880 int iSL3_hi = superlayer_boundaries[2];
881 int iSL5_lo = superlayer_boundaries[3]+1;
882 int iSL5_hi = superlayer_boundaries[3];
883
884 // Loop through seeds, looking for ones with only SL1 hits
885 for(unsigned int i=0; i<seeds.size(); i++){
886 if(!seeds[i].valid)continue;
887 vector<DCDCTrkHit*> &hits = seeds[i].hits;
888 if(hits.size()<1)continue;
889 bool has_SL1_hit = false;
890 bool has_SL3_hit = false;
891 bool has_SL5_hit = false;
892 for(unsigned int j=0; j<hits.size(); j++){
893 int ring = hits[j]->hit->wire->ring;
894 if(ring<iSL1_lo || ring>iSL1_hi)has_SL1_hit = true;
895 if(ring<iSL3_lo || ring>iSL3_hi)has_SL3_hit = true;
896 if(ring<iSL5_lo || ring>iSL5_hi)has_SL5_hit = true;
897 if(has_SL1_hit)break;
898 }
899 if(has_SL1_hit)continue;
900 if(has_SL3_hit && has_SL5_hit)continue;
901
902 // Seed contains hits only from superlayer 3 or only from superlayer 5
903 // Flag the seed as invalid so it is ignored later
904 seeds[i].valid = false;
905 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
905<<" "<<"Dropping seed "<<i<<": hits in oly SL3 or only SL5 (has_SL3_hit="<<has_SL3_hit<<", has_SL5_hit="<<has_SL5_hit<<")"<<endl;
906 }
907}
908
909//------------------
910// FilterCloneSeeds
911//------------------
912void DTrackCandidate_factory_CDC::FilterCloneSeeds(vector<DCDCSeed> &seeds)
913{
914 /// Look for clones of seeds and flag all but one as not valid.
915 /// Two tracks are considered clones if their trajectories are
916 /// really close.
917
918 /// The "closeness" of the two trajectories is determined using
919 /// two methods, either of which can flag the two seeds as clones.
920 ///
921 /// For the first method, we use the "asymmetry" of the circle
922 /// centers. The asymmetry is defined as the distance between
923 /// the two circle's centers divided by the sum of the magnitudes
924 /// of the circle's centers relative to the origin. If this is below
925 /// some threshold (say 0.05) then the tracks are seeds are
926 /// considered clones.
927 ///
928 /// For the second method, we check the distance between the
929 /// trajectories at a point near the last hit and a point near
930 /// one of the middle hits of one of the seeds. Since all
931 /// seeds pass through the origin, this effectively gives 3
932 /// points on the circle(s). If the total sum of the distances
933 /// between these points and the other seed's track are less than
934 /// some small value, the two are considered clones. Note that
935 /// currently, the choice of "middle" hits is done simply by
936 /// looking at the N/2th hit. This may not actually be near the
937 /// middle of the R values so a more sophisticated algorithm
938 /// may need to be implemented later.
939
940 for(unsigned int i=0; i<seeds.size(); i++){
941 DCDCSeed &seed1 = seeds[i];
942 if(!seed1.valid)continue;
943
944 // Calculate point on seed1 closest to last(furthest out) hit on seed1
945 DVector2 r1(seed1.fit.x0, seed1.fit.y0);
946 const DVector3 &wire_pos1 = seed1.hits[seed1.hits.size()-1]->hit->wire->origin;
947 DVector2 alpha1(wire_pos1.X(), wire_pos1.Y());
948 alpha1 -= r1;
949 alpha1 *= r1.Mod()/alpha1.Mod();
950
951 // Calculate point on seed1 at about the midpoint hit on seed1
952 const DVector3 &wire_pos2 = seed1.hits[seed1.hits.size()/2]->hit->wire->origin;
953 DVector2 alpha2(wire_pos2.X(), wire_pos2.Y());
954 alpha2 -= r1;
955 alpha2 *= r1.Mod()/alpha2.Mod();
956
957 for(unsigned int j=i+1; j<seeds.size(); j++){
958 DCDCSeed &seed2 = seeds[j];
959 if(!seed2.valid)continue;
960
961 // Flag to indicate these are clones
962 bool are_clones = false;
963
964 // If the two circle fits are significantly different (i.e.
965 // their x0, y0 coordinates aren't close), then skip assume
966 // these seeds aren't duplicates.
967 DVector2 r2(seed2.fit.x0, seed2.fit.y0);
968 DVector2 dr = r1-r2;
969 double asym = dr.Mod()/(r1.Mod()+r2.Mod());
970 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
970<<" "<<"asym="<<asym<<endl;
971 if(asym<0.05) are_clones = true;
972
973 if(!are_clones){
974 // Find total absolute distance from alpha1 and seed2 circle and
975 // alpha2 and seed2 circle.
976 DVector2 d1 = dr+alpha1;
977 DVector2 d2 = dr+alpha2;
978 double d = fabs(d1.Mod()-r2.Mod()) + fabs(d2.Mod() - r2.Mod());
979 if(d<1.5) are_clones = true;
980 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
980<<" "<<"d="<<d<<endl;
981 }
982 // Check that the charges are the same
983 if (are_clones){
984 if (seed1.fit.q!=seed2.fit.q) are_clones=false;
985
986 }
987 // Remove a clone in necessary
988 if(are_clones){
989 // These seeds appear to be clones of one another. Mark one as not valid.
990
991 if(seed1.fit.chisq<=seed2.fit.chisq){
992 seed2.valid = false;
993 }else{
994 seed1.valid = false;
995 }
996
997 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
997<<" "<<"Filtering clone circle seed (seed1.fit.chisq="<<seed1.fit.chisq<<" seed2.fit.chisq="<<seed2.fit.chisq<<endl;
998 }
999 }
1000 }
1001}
1002
1003//------------------
1004// AddStereoHits
1005//------------------
1006void DTrackCandidate_factory_CDC::AddStereoHits(vector<DCDCTrkHit*> &stereo_hits, DCDCSeed &seed)
1007{
1008 // To find the z coordinate, we look at the 2D projection of the
1009 // stereo wire and find the intersection point of that with the
1010 // circle found in FitSeed().
1011
1012 // Loop over stereo hits to find z-values for any that cross this circle
1013 for(unsigned int i=0; i<stereo_hits.size(); i++){
1014 DCDCTrkHit *trkhit = stereo_hits[i];
1015
1016 // Don't consider hits that are out of time with the seed
1017 if(fabs(trkhit->hit->tdrift-seed.tdrift_avg)>MAX_SEED_TIME_DIFF){
1018 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1018<<" "<<"Out of time stereo hit: seed.tdrift_avg="<<seed.tdrift_avg<<" tdrift="<<trkhit->hit->tdrift<<endl;
1019 continue;
1020 }
1021
1022 // Calculate intersection points between circle and stereo wire
1023 const DCDCWire *wire = trkhit->hit->wire;
1024 DVector2 r1(wire->origin.X(), wire->origin.Y());
1025 DVector2 r2(wire->udir.X(), wire->udir.Y());
1026 DVector2 R(seed.fit.x0, seed.fit.y0);
1027 double r2_mod2 = r2.Mod2();
1028 if(r2_mod2<1.0E-10){
1029 _DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1029<<" "<<"wire in CDC stereo hit list is not stereo!"<<endl;
1030 continue; // this must not be a stereo wire!
1031 }
1032
1033 double a = r2_mod2;
1034 double b = 2.0*r2*(r1-R);
1035 double c = r1.Mod2()-2.0*r1*R;
1036 double A = b*b - 4.0*a*c;
1037
1038 // Check that this wire intersects this circle
1039 if(A<0.0)continue; // line along wire does not intersect circle, ever.
1040
1041 // Calculate intersection points for two roots of alpha
1042 double B = sqrt(A);
1043 double alpha1 = (-b - B)/(2.0*a);
1044 double alpha2 = (-b + B)/(2.0*a);
1045
1046 if(DEBUG_LEVEL>15)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1046<<" "<<"alpha1="<<alpha1<<" alpha2="<<alpha2<<endl;
1047
1048 // At this point we must decide which value of alpha to use. The
1049 // proper way would likely involve either trying both possibilities
1050 // and taking the one that gave the better chi-sq for a line fit
1051 // of phi vs. z or looking at the surrounding axial layers
1052 // and using the value which puts the hit closest to those.
1053 // For now, we just use the value closest to zero (i.e. closest to
1054 // the center of the wire).
1055 //double alpha = sqrt(r2_mod2)*(fabs(alpha1)<fabs(alpha2) ? alpha1:alpha2);
1056
1057 // Now we must convert the alpha value into a z-value. To do this,
1058 // we use the known theta angle of the wire. The distance along the
1059 // wire from the center in 3D is given by:
1060 //
1061 // s = alpha/sin(theta_wire)
1062 //
1063 // The z coordinate of the hit is given by:
1064 //
1065 // z = z1 + s*z2
1066 //
1067 // where z1 is the z coordinate of the center of the wire and z2
1068 // is the z component of the direction of the wire. i.e.
1069 //
1070 // z2 = cos(theta_wire)
1071 //
1072 // This means sin(theta_wire) = sqrt(1 - (z2)^2)
1073 //double z2 = wire->udir.Z();
1074 //double s = alpha/fabs(sin(wire->stereo));
1075 // The above code unnecessarily multiples and divides by sin(theta_wire)
1076 double s=(fabs(alpha1)<fabs(alpha2) ? alpha1:alpha2);
1077
1078 //if(DEBUG_LEVEL>3)_DBG_<<"alpha="<<alpha<<" s="<<s<<" ring="<<wire->ring<<" straw="<<wire->straw<<" stereo="<<wire->stereo<<endl;
1079 if(DEBUG_LEVEL>15)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1079<<" "<<"s="<<s<<" ring="<<wire->ring<<" straw="<<wire->straw<<" stereo="<<wire->stereo<<endl;
1080 if(fabs(s) > wire->L/2.0)continue; // if wire doesn't cross circle, skip hit
1081
1082
1083 // Determine 3d space point for this hit
1084 DCDCTrkHit mytrkhit = *trkhit; // we need to make a copy since x_stereo,... are unique based on the cadidate
1085 DVector3 pos = wire->origin + s*wire->udir;
1086 mytrkhit.x_stereo = pos.X();
1087 mytrkhit.y_stereo = pos.Y();
1088 mytrkhit.z_stereo = pos.Z();
1089
1090 // Verify that this hit is reasonably close to the axial hits in X and Y
1091 // If not, then drop this hit for this seed.
1092 double phi_diff = fabs(pos.Phi() - seed.phi_avg);
1093 if(phi_diff>M_PI3.14159265358979323846)phi_diff = 2.0*M_PI3.14159265358979323846 - phi_diff;
1094 phi_diff*=57.3; // convert to degrees
1095 if(fabs(phi_diff)>MAX_SEED_LINK_ANGLE){ // yes, the fabs is needed ...
1096 if(DEBUG_LEVEL>4)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1096<<" "<<"rejecting stereo hit at phi="<<pos.Phi()<<" for being too far away from axial hits(phi_diff="<<phi_diff<<")"<<endl;
1097 continue;
1098 }
1099 mytrkhit.phi_stereo = atan2(mytrkhit.y_stereo-R.Y(), mytrkhit.x_stereo-R.X());
1100 R*=-1.0; // make R point from center of circle to beamline instead of other way around
1101 if(DEBUG_LEVEL>15){
1102 _DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1102<<" "<<" --- wire->udir X, Y, Z = "<<wire->udir.X()<<", "<<wire->udir.Y()<<", "<<wire->udir.Z()<<endl;
1103 _DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1103<<" "<<" -- ring="<<wire->ring<<" trkhit->z_stereo="<<trkhit->z_stereo<<" trkhit->y_stereo="<<trkhit->y_stereo<<" trkhit->x_stereo="<<trkhit->x_stereo<<endl;
1104 _DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1104<<" "<<" -- phi_stereo="<<trkhit->phi_stereo<<" R.Phi()="<<R.Phi()<<" (X,Y)=("<<R.X()<<", "<<R.Y()<<")"<<endl;
1105 _DBG__std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1105<<std::endl;
1106 }
1107 mytrkhit.phi_stereo -= R.Phi(); // make angle relative to beamline
1108
1109 // We want this to go either from 0 to +2pi for positive charge, or
1110 // 0 to -2pi for negative.
1111 double phi_hi = seed.fit.q>0.0 ? +2.0*M_PI3.14159265358979323846:0.0;
1112 double phi_lo = seed.fit.q>0.0 ? 0.0:-2.0*M_PI3.14159265358979323846;
1113 while(mytrkhit.phi_stereo<phi_lo){
1114 mytrkhit.phi_stereo+=2.0*M_PI3.14159265358979323846;
1115 }
1116 while(mytrkhit.phi_stereo>phi_hi){
1117 mytrkhit.phi_stereo-=2.0*M_PI3.14159265358979323846;
1118 }
1119 mytrkhit.flags |= VALID_STEREO;
1120 seed.stereo_hits.push_back(mytrkhit);
1121 if(DEBUG_LEVEL>10)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1121<<" "<<"Adding CDC stereo hit: ring="<<mytrkhit.hit->wire->ring<<" straw="<<mytrkhit.hit->wire->straw<<endl;
1122 }
1123 if(DEBUG_LEVEL>5)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1123<<" "<<"Num stereo hits: "<<seed.stereo_hits.size()<<endl;
1124}
1125
1126//------------------
1127// FindThetaZ
1128//------------------
1129void DTrackCandidate_factory_CDC::FindThetaZ(DCDCSeed &seed)
1130{
1131 // Decide whether to use the histogramming method or the
1132 // straight track method base on the transverse momentum
1133if(DEBUG_LEVEL>2)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1133<<" "<<"p_trans:"<<seed.fit.p_trans<<endl;
1134//if(seed.fit.p_trans<3.0)
1135 if (true)
1136 {
1137 FindTheta(seed, TARGET_Z_MIN, TARGET_Z_MAX);
1138 FindZ(seed, seed.theta_min, seed.theta_max);
1139 }else{
1140 FindThetaZStraightTrack(seed);
1141
1142 // The value of p_trans that was originally determined
1143 // by DQuickFit::FitCircleStraightTrack() can often lead
1144 // to a total momentum larger than the beam momentum.
1145 // Re-estimate p_trans limiting the search with the
1146 // theta we just found.
1147 seed.fit.SearchPtrans(9.0*sin(seed.theta), 0.1);
1148 //seed.fit.QuickPtrans();
1149 if(DEBUG_LEVEL>2)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1149<<" "<<"p_trans from search:"<<seed.fit.p_trans<<" (ptot="<<seed.fit.p_trans/sin(seed.theta)<<")"<<endl;
1150 }
1151
1152 // If z_vertex is not inside the target limits, then flag this
1153 // seed as invalid.
1154 /*
1155 if(seed.z_vertex<TARGET_Z_MIN || seed.z_vertex>TARGET_Z_MAX){
1156 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;
1157 seed.valid=false;
1158 }
1159
1160 */
1161 return;
1162}
1163
1164//------------------
1165// FindThetaZStraightTrack
1166//------------------
1167void DTrackCandidate_factory_CDC::FindThetaZStraightTrack(DCDCSeed &seed)
1168{
1169 /// In the case of high momentum tracks, the values of phi_stereo will
1170 /// all be close together (and near 0 or +-2pi). For these cases,
1171 /// we calculate theta from r and z points which should be pretty linear.
1172
1173 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1173<<" "<<"Fitting theta and z for high transverse momentum track (p_trans="<<seed.fit.p_trans<<")"<<endl;
1
Taking false branch
→
1174
1175 // Do a linear regression of R vs. Z for stereo hits. First, we make
1176 // a list of the hits we'll use.
1177 vector<double> r;
1178 vector<double> z;
1179 for(unsigned int i=0; i<seed.stereo_hits.size(); i++){
2
←
Loop condition is false. Execution continues on line 1191
→
1180 DCDCTrkHit *trkhit = &seed.stereo_hits[i];
1181 if(!trkhit->flags&VALID_STEREO)continue;
1182
1183 //double R = sqrt(pow(trkhit->x_stereo,2.0) + pow(trkhit->y_stereo,2.0));
1184 double R=sqrt(trkhit->x_stereo*trkhit->x_stereo
1185 +trkhit->y_stereo*trkhit->y_stereo);
1186 r.push_back(R);
1187 z.push_back(trkhit->z_stereo);
1188 }
1189
1190 // Add center of target as a point to constrain the fit a little more
1191 r.push_back(0.0);
1192 z.push_back(TARGET_Z);
1193
1194 // Calculate average z and r
1195 double Ravg=0.0, Zavg=0.0;
1196 for(unsigned int i=0; i<r.size(); i++){
3
←
Loop condition is true. Entering loop body
→
1197 Ravg += r[i];
4
←
Dereference of null pointer
1198 Zavg += z[i];
1199 }
1200 Ravg /= (double)r.size();
1201 Zavg /= (double)z.size();
1202
1203 // Now, do the regression
1204 double Szr=0.0, Szz=0.0, Srr=0.0;;
1205 for(unsigned int i=0; i<r.size(); i++){
1206 if(DEBUG_LEVEL>4)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1206<<" "<<"r="<<r[i]<<"\t z="<<z[i]<<" Szz_i="<<pow((z[i] - Zavg), 2.0)<<" Srr_i="<<pow((r[i] - Ravg),2.0)<<endl;
1207 double dz=z[i] - Zavg;
1208 double dr=r[i] - Ravg;
1209 /*
1210 Szr += (z[i] - Zavg)*(r[i] - Ravg);
1211 Szz += pow((z[i] - Zavg), 2);
1212 Srr += pow((r[i] - Ravg), 2);
1213 */
1214 Szr += dz*dr;
1215 Szz += dz*dz;
1216 Srr += dr*dr;
1217 }
1218
1219 if(Szz>Srr){
1220 // track is more horizontal
1221 seed.z_vertex = Zavg - Ravg*Szz/Szr;
1222 seed.theta = atan2(Szr, Szz);
1223 }else{
1224 // track is more vertical
1225 seed.z_vertex = Zavg - Ravg*Szr/Srr;
1226 seed.theta = M_PI_21.57079632679489661923 - atan2(Szr, Srr);
1227 }
1228 if(seed.theta<0.0){
1229 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1229<<" "<<"theta less than 0 ("<<seed.theta<<") this simply shoudn't happen!"<<endl;
1230 seed.theta+=M_PI3.14159265358979323846;
1231 }
1232 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1232<<" "<<"theta="<<seed.theta<<" z_vertex="<<seed.z_vertex<<endl;
1233 if(DEBUG_LEVEL>4){
1234 _DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1234<<" "<<"Szz="<<Szz<<" Srr="<<Srr<<endl;
1235 double z_vertex = Zavg - Ravg*Szz/Szr;
1236 double theta = atan2(Szr, Szz);
1237 _DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1237<<" "<<"for Szz>Srr : theta="<<theta<<" z_vertex="<<z_vertex<<endl;
1238 z_vertex = Zavg - Ravg*Szr/Srr;
1239 theta = M_PI_21.57079632679489661923 - atan2(Szr, Srr);
1240 _DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1240<<" "<<"for Szz<=Srr : theta="<<theta<<" z_vertex="<<z_vertex<<endl;
1241 }
1242}
1243
1244//------------------
1245// FindTheta
1246//------------------
1247void DTrackCandidate_factory_CDC::FindTheta(DCDCSeed &seed, double target_z_min, double target_z_max)
1248{
1249 /// Find the theta value using the valid stereo hits from <i>seed</i>. The values
1250 /// for phi_stereo and z_stereo are assumed to be valid as is the status of the
1251 /// VALID_STEREO bit in flags. The value of seed.fit.r0 is also used to calculate
1252 /// theta.
1253 ///
1254 /// This uses a histogramming technique that looks at the overlaps of the
1255 /// angle ranges subtended by each hit between the given target limits.
1256 /// The overlaps usually lead to a range of values for theta. The limits
1257 /// of these are stored in the theta_min and theta_max fields of the seed.
1258 /// The centroid of the range is stored in the theta field.
1259
1260 // We use a simple array to store our histogram here. We don't want to use
1261 // ROOT histograms because they are not thread safe.
1262 unsigned int Nbins = 1000;
1263 unsigned int hist[Nbins];
1264 for(unsigned int i=0; i<Nbins; i++)hist[i] = 0; // clear histogram
1265 double bin_width = 2.0*M_PI3.14159265358979323846/(double)Nbins;
1266 double hist_low_limit = -M_PI3.14159265358979323846; // lower edge of histogram limits
1267
1268 // Loop over CDC hits, filling the histogram
1269 double r0 = seed.fit.r0;
1270 for(unsigned int i=0; i<seed.stereo_hits.size(); i++){
1271 DCDCTrkHit *trkhit = &seed.stereo_hits[i];
1272 if(!trkhit->flags&VALID_STEREO)continue;
1273
1274 // Calculate upper and lower limits in theta
1275 double alpha = r0*trkhit->phi_stereo;
1276 if(seed.fit.q<0.0)alpha = -alpha;
1277 double tmin = atan2(alpha, trkhit->z_stereo - target_z_min);
1278 double tmax = atan2(alpha, trkhit->z_stereo - target_z_max);
1279 if(tmin>tmax){
1280 double tmp = tmin;
1281 tmin=tmax;
1282 tmax=tmp;
1283 }
1284 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1284<<" "<<" -- phi_stereo="<<trkhit->phi_stereo<<" z_stereo="<<trkhit->z_stereo<<" alpha="<<alpha<<endl;
1285 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1285<<" "<<" -- tmin="<<tmin<<" tmax="<<tmax<<endl;
1286
1287 // Find index of bins corresponding to tmin and tmax
1288 unsigned int imin = (unsigned int)floor((tmin-hist_low_limit)/bin_width);
1289 unsigned int imax = (unsigned int)floor((tmax-hist_low_limit)/bin_width);
1290
1291 // If entire range of this hit is outside of the histogram limit
1292 // then discard this hit.
1293 if(imax<0 || imin>=Nbins)continue;
1294
1295 // Clip limits of bin range to our histogram limits
1296 if(imin<0)imin=0;
1297 if(imin>=Nbins)imin=Nbins-1;
1298 if(imax<0)imax=0;
1299 if(imax>=Nbins)imax=Nbins-1;
1300
1301 // Increment all bins between imin and imax
1302 for(unsigned int j=imin; j<=imax; j++)hist[j]++;
1303 }
1304
1305 // Look for the indexes of the plateau
1306 unsigned int istart=0;
1307 unsigned int iend=0;
1308 for(unsigned int i=1; i<Nbins; i++){
1309 if(hist[i]>hist[istart]){
1310 istart = i;
1311 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1311<<" "<<" -- istart="<<istart<<" (theta="<<hist_low_limit + bin_width*(0.5+(double)istart)<<" , N="<<hist[i]<<")"<<endl;
1312 }
1313 if(hist[i] == hist[istart])iend = i;
1314 }
1315
1316 // If there are no entries in the histogram, then flag this seed as invalid
1317 if(hist[istart]==0.0)seed.valid=false;
1318
1319 // Calculate theta limits
1320 seed.theta_min = hist_low_limit + bin_width*(0.5+(double)istart);
1321 seed.theta_max = hist_low_limit + bin_width*(0.5+(double)iend);
1322 seed.theta = (seed.theta_max + seed.theta_min)/2.0;
1323 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1323<<" "<<"istart="<<istart<<" iend="<<iend<<" theta_min="<<seed.theta_min<<" theta_max="<<seed.theta_max<<endl;
1324}
1325
1326//------------------
1327// FindZ
1328//------------------
1329void DTrackCandidate_factory_CDC::FindZ(DCDCSeed &seed, double theta_min, double theta_max)
1330{
1331 /// Find the z value of the vertex using the valid stereo hits from <i>seed</i>. The values
1332 /// for phi_stereo and z_stereo are assumed to be valid as is the status of the
1333 /// VALID_STEREO bit in flags.
1334 ///
1335 /// This uses a histogramming technique that looks at the overlaps of the
1336 /// z ranges subtended by each hit between the given theta limits.
1337 /// The overlaps usually lead to a range of values for z_vertex. The limits
1338 /// of these are stored in the z_min and z_max fields of the seed.
1339 /// The centroid of the range is stored in the z_vertex field.
1340
1341 // We use a simple array to store our histogram here. We don't want to use
1342 // ROOT histograms because they are not thread safe.
1343 unsigned int Nbins = 300;
1344 unsigned int hist[Nbins];
1345 for(unsigned int i=0; i<Nbins; i++)hist[i] = 0; // clear histogram
1346 double bin_width = 0.5; // bins are 5mm
1347 double hist_low_limit = 0.0; // lower edge of histogram limits
1348
1349 // Loop over CDC hits, filling the histogram
1350 double r0 = seed.fit.r0;
1351 double tan_alpha_min = tan(theta_min)/r0;
1352 double tan_alpha_max = tan(theta_max)/r0;
1353 for(unsigned int i=0; i<seed.stereo_hits.size(); i++){
1354 DCDCTrkHit *trkhit = &seed.stereo_hits[i];
1355 if(!trkhit->flags&VALID_STEREO)continue;
1356
1357 // Calculate upper and lower limits in z
1358 double q_sign = seed.fit.q>0.0 ? +1.0:-1.0;
1359 double zmin = trkhit->z_stereo - q_sign*trkhit->phi_stereo/tan_alpha_min;
1360 double zmax = trkhit->z_stereo - q_sign*trkhit->phi_stereo/tan_alpha_max;
1361 if(zmin>zmax){
1362 double tmp = zmin;
1363 zmin=zmax;
1364 zmax=tmp;
1365 }
1366 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1366<<" "<<" -- phi_stereo="<<trkhit->phi_stereo<<" z_stereo="<<trkhit->z_stereo<<endl;
1367 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1367<<" "<<" -- zmin="<<zmin<<" zmax="<<zmax<<endl;
1368
1369 // Find index of bins corresponding to tmin and tmax
1370 unsigned int imin = (unsigned int)floor((zmin-hist_low_limit)/bin_width);
1371 unsigned int imax = (unsigned int)floor((zmax-hist_low_limit)/bin_width);
1372
1373 // If entire range of this hit is outside of the histogram limit
1374 // then discard this hit.
1375 if(imax<=0 || imin>=Nbins)continue;
1376
1377 // Clip limits of bin range to our histogram limits
1378 if(imin<0)imin=0;
1379 if(imin>=Nbins)imin=Nbins-1;
1380 if(imax<0)imax=0;
1381 if(imax>=Nbins)imax=Nbins-1;
1382
1383 // Increment all bins between imin and imax
1384 for(unsigned int j=imin; j<=imax; j++)hist[j]++;
1385 }
1386
1387 // Look for the indexes of the plateau
1388 unsigned int istart=0;
1389 unsigned int iend=0;
1390 for(unsigned int i=1; i<Nbins; i++){
1391 if(hist[i]>hist[istart]){
1392 istart = i;
1393 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1393<<" "<<" -- istart="<<istart<<" (z="<<hist_low_limit + bin_width*(0.5+(double)istart)<<" , N="<<hist[i]<<")"<<endl;
1394 }
1395 if(hist[i] == hist[istart])iend = i;
1396 }
1397
1398 // If there are no entries in the histogram, then flag this seed as invalid
1399 if(hist[istart]==0.0)seed.valid=false;
1400
1401 // Calculate z limits
1402 seed.z_min = hist_low_limit + bin_width*(0.5+(double)istart);
1403 seed.z_max = hist_low_limit + bin_width*(0.5+(double)iend);
1404 seed.z_vertex = (seed.z_max + seed.z_min)/2.0;
1405 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1405<<" "<<"istart="<<istart<<" iend="<<iend<<" z_min="<<seed.z_min<<" z_max="<<seed.z_max<<" hits[istart]="<<hist[istart]<<endl;
1406}
1407
1408//------------------
1409// NumEligibleSeedHits
1410//------------------
1411int DTrackCandidate_factory_CDC::NumEligibleSeedHits(vector<DCDCTrkHit*> &hits)
1412{
1413 int N=0;
1414 for(unsigned int i=0; i<hits.size(); i++){
1415 if((hits[i]->flags & (USED | CANT_BE_IN_SEED | OUT_OF_TIME)) == 0)N++;
1416 }
1417
1418 return N;
1419}
1420
1421//------------------
1422// DCDCSeed::DCDCSeed
1423//------------------
1424DTrackCandidate_factory_CDC::DCDCSeed::DCDCSeed()
1425{
1426 phi_avg = 0.0;
1427 tdrift_avg = 0.0;
1428 linked = false;
1429 valid = false;
1430 theta = 0.0;
1431 z_vertex = 0.0;
1432 q = 0.0;
1433 theta_min = 0.0;
1434 theta_max = 0.0;
1435 z_min = 0.0;
1436 z_max = 0.0;
1437}
1438
1439//------------------
1440// DCDCSeed::Merge
1441//------------------
1442void DTrackCandidate_factory_CDC::DCDCSeed::Merge(DCDCSeed& seed)
1443{
1444 // Need avg. phi based on hits from both seeds
1445 double x = (double)hits.size()*cos(phi_avg);
1446 double y = (double)hits.size()*sin(phi_avg);
1447
1448 for(unsigned int i=0; i<seed.hits.size(); i++){
1449 hits.push_back(seed.hits[i]);
1450
1451 x += cos(seed.hits[i]->hit->wire->phi);
1452 y += sin(seed.hits[i]->hit->wire->phi);
1453 }
1454
1455 phi_avg = atan2(y,x);
1456 if(phi_avg<0.0)phi_avg += 2.0*M_PI3.14159265358979323846;
1457}
1458
1459//------------------
1460// DCDCSeed::MinDist2
1461//------------------
1462double DTrackCandidate_factory_CDC::DCDCSeed::MinDist2(DCDCSeed& seed)
1463{
1464 /// Returns the minimum distance squared between this and the given seed.
1465 /// Only the first and last hits of each seed's hit list are used
1466 /// to calculate a maximum of 4 distances (minimum of 1) of which, the
1467 /// smallest is returned.
1468 if(hits.size()==0 || seed.hits.size()==0){
1469 _DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1469<<" "<<"Number of seed hits 0! (Nthis="<<hits.size()<<" ,Nseed="<<seed.hits.size()<<")"<<endl;
1470 return 1.0E10;
1471 }
1472
1473 // This is kind of messy looking, but it avoids extra calls to
1474 // DCDCTrkHit::Dist2 that aren't necessary.
1475 double d2min = hits[0]->Dist2(seed.hits[0]);
1476 if(hits.size()>0){
1477 double d2 = hits[hits.size()-1]->Dist2(seed.hits[0]);
1478 if(d2<d2min)d2min = d2;
1479
1480 if(seed.hits.size()>0){
1481 double d2 = hits[hits.size()-1]->Dist2(seed.hits[seed.hits.size()-1]);
1482 if(d2<d2min)d2min = d2;
1483 }
1484 }
1485 if(seed.hits.size()>0){
1486 double d2 = hits[0]->Dist2(seed.hits[seed.hits.size()-1]);
1487 if(d2<d2min)d2min = d2;
1488 }
1489
1490 return d2min;
1491}
1492
1493//------------------
1494// DCDCSeed::FindAverageBz
1495//------------------
1496double DTrackCandidate_factory_CDC::DCDCSeed::FindAverageBz( const DMagneticFieldMap *bfield
1497)
1498{
1499 //return 2.0;
1500 if(!bfield)return 0.0;
1501
1502 double Bz_sum=0.0;
1503 for(unsigned int i=0; i<stereo_hits.size(); i++){
1504 DCDCTrkHit *hit = &stereo_hits[i];
1505 Bz_sum += bfield->GetBz(hit->x_stereo,hit->y_stereo,
1506 hit->z_stereo);
1507 }
1508
1509 return Bz_sum/(double)(stereo_hits.size());
1510}
1511
1512
1513//------------------
1514// FindThetaZRegression
1515//------------------
1516// Linear regression to find tan(lambda) and z_vertex
1517jerror_t DTrackCandidate_factory_CDC::FindThetaZRegression(DCDCSeed &seed){
1518
1519 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1519<<" "<<"Finding theta and z via linear regression method."<<endl;
1520
1521 if (seed.fit.normal.Mag()==0.) return VALUE_OUT_OF_RANGE;
1522 // Vector of intersections between the circles of the measurements and the plane intersecting the Riemann surface
1523 vector<DVector3_with_perp>intersections;
1524 DVector3_with_perp beam(0,0,TARGET_Z);
1525 intersections.push_back(beam);
1526
1527 // CDC stereo hits
1528 for (unsigned int m=0;m<seed.stereo_hits.size();m++){
1529 unsigned int numbits=8*sizeof(unsigned int);
1530 seed.HitBitPattern[seed.stereo_hits[m].index/numbits]
1531 |=1<<seed.stereo_hits[m].index%numbits;
1532
1533 DCDCTrkHit *trkhit=&seed.stereo_hits[m];
1534 double R2=trkhit->x_stereo*trkhit->x_stereo+trkhit->y_stereo*trkhit->y_stereo;
1535
1536 DVector3_with_perp intersection;
1537 DVector3 N=seed.fit.normal;
1538 //double c0=seed.fit.c_origin;
1539 double A=seed.fit.c_origin+R2*N.Z();
1540 double B=N.Perp();
1541 double C=B*R2-A*A;
1542
1543 if (C>=0) {
1544 double sqrtC=sqrt(C);
1545 double x1=(-N.X()*A+N.Y()*sqrtC)/B;
1546 double y1=(-N.Y()*A-N.X()*sqrtC)/B;
1547 double x2=(-N.X()*A-N.Y()*sqrtC)/B;
1548 double y2=(-N.Y()*A+N.X()*sqrtC)/B;
1549
1550 if (fabs(trkhit->y_stereo-y1)<fabs(trkhit->y_stereo-y2)){
1551 intersection.SetX(x1);
1552 intersection.SetY(y1);
1553 }
1554 else{
1555 intersection.SetX(x2);
1556 intersection.SetY(y2);
1557 }
1558 intersection.SetZ(trkhit->z_stereo);
1559
1560 intersections.push_back(intersection);
1561 if(DEBUG_LEVEL>5)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1561<<" "<<"Adding CDC hit "<<m<<" z="<<intersection.z()<<endl;
1562
1563 }
1564
1565 }
1566
1567 // The SortIntersections method requires the perp field to be valid
1568 // for all objects. Go through and fill in those fields now.
1569 for(unsigned int i=0; i<intersections.size(); i++)intersections[i].CalcPerp();
1570
1571 // Now, sort the entries
1572 sort(intersections.begin(),intersections.end(),SortIntersections);
1573
1574 // Compute the arc lengths between (0,0) and (xi,yi)
1575 vector<double>arclengths(intersections.size());
1576 vector<double>var_z(intersections.size());
1577 arclengths[0]=0.;
1578 //double temp=1.6/sin(6.*M_PI/180.);
1579 //double varz=temp*temp/12.;
1580 double two_rc=2.*seed.fit.r0;
1581 for (unsigned int m=1;m<arclengths.size();m++){
1582 double diffx=intersections[m].x()-intersections[0].x();
1583 double diffy=intersections[m].y()-intersections[0].y();
1584 double chord=sqrt(diffx*diffx+diffy*diffy);
1585 double ratio=chord/two_rc;
1586 double s=(ratio<1.)?two_rc*asin(ratio):M_PI_21.57079632679489661923*two_rc;
1587
1588 arclengths[m]=s;
1589 //var_z[m]=varz;
1590 var_z[m]=19.52;
1591 }
1592
1593 //Linear regression to find z0, tanl
1594 double tanl=0.,z0=0.;
1595 if (arclengths.size()>2){ // Do fit only if have more than one measurement
1596 //var_z[0]=30.*30./12.; // Use length of the target
1597 var_z[0]=75.;
1598 double sumv=0.,sumx=0.;
1599 double sumy=0.,sumxx=0.,sumxy=0.;
1600 for (unsigned int m=0;m<intersections.size();m++){
1601 double temp=1./var_z[m];
1602 sumv+=temp;
1603 sumx+=arclengths[m]*temp;
1604 sumy+=intersections[m].z()*temp;;
1605 sumxx+=arclengths[m]*arclengths[m]*temp;
1606 sumxy+=arclengths[m]*intersections[m].z()*temp;
1607 }
1608 double Delta=sumv*sumxx-sumx*sumx;
1609 if (Delta==0.) return VALUE_OUT_OF_RANGE;
1610
1611 tanl=(sumv*sumxy-sumx*sumy)/Delta;
1612 z0=(sumxx*sumy-sumx*sumxy)/Delta;
1613 }
1614 else if(arclengths.size()==2){
1615 z0=TARGET_Z;
1616 tanl=(intersections[1].z()-z0)/arclengths[1];
1617 }else{
1618 if(DEBUG_LEVEL>5)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1618<<" "<<"Fit failed for theta-z via regressionz due to too few hits with z-info"<<endl;
1619 return VALUE_OUT_OF_RANGE;
1620 }
1621
1622 if (z0>TARGET_Z_MAX || z0<TARGET_Z_MIN){
1623 if(DEBUG_LEVEL>5)_DBG_std::cerr<<"DTrackCandidate_factory_CDC.cc"<<":"<<
1623<<" "<<"Fit failed for theta-z via regressionz value out of target range (z="<<z0<<")"<<endl;
1624 return VALUE_OUT_OF_RANGE;
1625 }
1626
1627 seed.theta=M_PI_21.57079632679489661923-atan(tanl);
1628 seed.z_vertex=z0;
1629
1630 return NOERROR;
1631}
1632
1633// Get the position and momentum at a fixed radius from the beam line
1634jerror_t DTrackCandidate_factory_CDC::GetPositionAndMomentum(DCDCSeed &seed,
1635 DVector3 &pos,
1636 DVector3 &mom){
1637 // Direction tangent and transverse momentum
1638 double tanl=tan(M_PI_21.57079632679489661923-seed.theta);
1639 double pt=seed.fit.p_trans*fabs(seed.FindAverageBz(bfield))/2.0;
1640
1641 // Squared radius of cylinder outside start counter but inside CDC inner
1642 // radius
1643 double r2=90.0;
1644
1645 // Circle parameters
1646 double xc=seed.fit.x0;
1647 double yc=seed.fit.y0;
1648 double rc=seed.fit.r0;
1649 double rc2=rc*rc;
1650 double xc2=xc*xc;
1651 double yc2=yc*yc;
1652 double xc2_plus_yc2=xc2+yc2;
1653
1654 // variables needed for intersection of circles
1655 double a=(r2-xc2_plus_yc2-rc2)/(2.*rc);
1656 double temp1=yc*sqrt(xc2_plus_yc2-a*a);
1657 double temp2=xc*a;
1658 double cosphi_plus=(temp2+temp1)/xc2_plus_yc2;
1659 double cosphi_minus=(temp2-temp1)/xc2_plus_yc2;
1660
1661 // Check for intersections
1662 if(!isfinite(temp1) || !isfinite(temp2)){
1663 // We did not find an intersection between the two circles, so return
1664 // sensible defaults for pos and
1665 double my_phi=seed.fit.phi;
1666 double sinphi=sin(my_phi);
1667 double cosphi=cos(my_phi);
1668 double D=-seed.fit.q*rc-xc/sinphi;
1669 pos.SetXYZ(-D*sinphi,D*cosphi,seed.z_vertex);
1670 mom.SetXYZ(pt*cosphi,pt*sinphi,pt*tanl);
1671
1672 return NOERROR;
1673 }
1674
1675 // if we have intersections, find both solutions
1676 double phi_plus=-acos(cosphi_plus);
1677 double phi_minus=-acos(cosphi_minus);
1678 double x_plus=xc+rc*cosphi_plus;
1679 double x_minus=xc+rc*cosphi_minus;
1680 double y_plus=yc+rc*sin(phi_plus);
1681 double y_minus=yc+rc*sin(phi_minus);
1682
1683 // if the resulting radial position on the circle from the fit does not agree
1684 // with the radius to which we are matching, we have the wrong sign for phi+
1685 // or phi-
1686 double r2_plus=x_plus*x_plus+y_plus*y_plus;
1687 double r2_minus=x_minus*x_minus+y_minus*y_minus;
1688 if (fabs(r2-r2_plus)>EPS1e-3){
1689 phi_plus*=-1.;
1690 y_plus=yc+rc*sin(phi_plus);
1691 }
1692 if (fabs(r2-r2_minus)>EPS1e-3){
1693 phi_minus*=-1.;
1694 y_minus=yc+rc*sin(phi_minus);
1695 }
1696
1697 // Choose phi- or phi+ depending on proximity to one of the cdc hits
1698 double xwire=seed.hits[0]->hit->wire->origin.x();
1699 double ywire=seed.hits[0]->hit->wire->origin.y();
1700 double dx=x_minus-xwire;
1701 double dy=y_minus-ywire;
1702 double d2_minus=dx*dx+dy*dy;
1703 dx=x_plus-xwire;
1704 dy=y_plus-ywire;
1705 double d2_plus=dx*dx+dy*dy;
1706 if (d2_plus>d2_minus){
1707 phi_minus*=-1.;
1708 if (seed.fit.q<0) phi_minus+=M_PI3.14159265358979323846;
1709 double dphi=M_PI_21.57079632679489661923-phi_minus-seed.fit.phi;
1710 while (dphi>2.*M_PI3.14159265358979323846) dphi-=2*M_PI3.14159265358979323846;
1711 while (dphi<-2.*M_PI3.14159265358979323846) dphi+=2*M_PI3.14159265358979323846;
1712 if (dphi<-M_PI3.14159265358979323846) dphi+=2*M_PI3.14159265358979323846;
1713 if (dphi>M_PI3.14159265358979323846) dphi-=2*M_PI3.14159265358979323846;
1714 pos.SetXYZ(x_minus,y_minus,seed.z_vertex+seed.fit.q*rc*dphi*tanl);
1715 mom.SetXYZ(pt*sin(phi_minus),pt*cos(phi_minus),pt*tanl);
1716
1717 }
1718 else{
1719 phi_plus*=-1.;
1720 if (seed.fit.q<0) phi_plus+=M_PI3.14159265358979323846;
1721 double dphi=M_PI_21.57079632679489661923-phi_plus-seed.fit.phi;
1722 while (dphi>2.*M_PI3.14159265358979323846) dphi-=2*M_PI3.14159265358979323846;
1723 while (dphi<-2.*M_PI3.14159265358979323846) dphi+=2*M_PI3.14159265358979323846;
1724 if (dphi<-M_PI3.14159265358979323846) dphi+=2*M_PI3.14159265358979323846;
1725 if (dphi>M_PI3.14159265358979323846) dphi-=2*M_PI3.14159265358979323846;
1726 pos.SetXYZ(x_plus,y_plus,seed.z_vertex+seed.fit.q*rc*dphi*tanl);
1727 mom.SetXYZ(pt*sin(phi_plus),pt*cos(phi_plus),pt*tanl);
1728
1729 }
1730
1731 return NOERROR;
1732}