File: | libraries/TRACKING/DTrackCandidate_factory_CDC.cc |
Location: | line 1207, column 13 |
Description: | Dereference of null pointer |
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> | |||
9 | using 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 | ||||
26 | inline 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 | ||||
38 | class 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 | ||||
50 | inline 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 | ||||
55 | inline 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 | } | |||
62 | inline 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 | } | |||
66 | inline 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 | } | |||
70 | inline 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 | //------------------ | |||
78 | jerror_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 | //------------------ | |||
114 | jerror_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 | //------------------ | |||
150 | jerror_t DTrackCandidate_factory_CDC::erun(void) | |||
151 | { | |||
152 | return NOERROR; | |||
153 | } | |||
154 | ||||
155 | //------------------ | |||
156 | // fini | |||
157 | //------------------ | |||
158 | jerror_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 | //------------------ | |||
167 | jerror_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 | //------------------ | |||
366 | jerror_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 | //------------------ | |||
453 | void 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 | //------------------ | |||
533 | void 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 | //------------------ | |||
602 | void 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 | //------------------ | |||
737 | bool 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 | //------------------ | |||
818 | void 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 | //------------------ | |||
873 | void 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 | //------------------ | |||
912 | void 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 | //------------------ | |||
1006 | void 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 | //------------------ | |||
1129 | void 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 | |||
1133 | if(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 | //------------------ | |||
1167 | void 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; | |||
| ||||
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++){ | |||
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++){ | |||
1197 | Ravg += r[i]; | |||
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 | //------------------ | |||
1247 | void 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 | //------------------ | |||
1329 | void 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 | //------------------ | |||
1411 | int 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 | //------------------ | |||
1424 | DTrackCandidate_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 | //------------------ | |||
1442 | void 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 | //------------------ | |||
1462 | double 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 | //------------------ | |||
1496 | double 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 | |||
1517 | jerror_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 | |||
1634 | jerror_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 | } |