CGEM BOSS 6.6.5.h
BESIII Offline Software System
Loading...
Searching...
No Matches
FTSuperLayer.cxx
Go to the documentation of this file.
1// -*- C++ -*-
2//
3// Package: MdcFastTrkAlg
4// Module: FTSuperLayer
5//
6// Description: super-layer class for MdcFastTrkAlg
7//
8#include "GaudiKernel/StatusCode.h"
9#include "GaudiKernel/Bootstrap.h"
10
11#include "MdcGeomSvc/IMdcGeomSvc.h"
12#include "MdcGeomSvc/MdcGeoLayer.h"
13
20#include <iostream>
21
23
24//...Globals...
25/*const unsigned int
26FTSuperLayer::_neighborsMask[6] = {452,420,340,172,92,60};*/
27// [0] = 452 = FTWireHitAppended + FTWireNeighbor345
28// [1] = 420 = FTWireHitAppended + FTWireNeighbor245
29// definitions of [2] = 340 = FTWireHitAppended + FTWireNeighbor135
30// _neighborsMask[] [3] = 172 = FTWireHitAppended + FTWireNeighbor024
31// [4] = 92 = FTWireHitAppended + FTWireNeighbor013
32// [5] = 60 = FTWireHitAppended + FTWireNeighbor012
33const float FTSuperLayer::_maxDphi[11] =
34 {11.9, 9.72, 6.26, 4.86, 3.81, 2.37, 2.08, 1.76, 1.53, 2.33, 1.88};
35
36void
38{
39 if (_wireHits.length()){
40 register FTWire ** hptr = _wireHits.firstPtr();
41 FTWire ** const last = _wireHits.lastPtr();
42 do{(**hptr).state(FTWireHitInvalid);}while((hptr++)!=last);
43 _wireHits.clear();
44 }
45 _singleHits.clear();
46 if (_segments.length()) _segments.deleteAll();
47 if (_complecated_segments) _complecated_segments->deleteAll();
48}
49
50
51void
53{
54 clustering();
55 FTList<FTSegment *> inner_short(10);
56 FTList<FTSegment *> outer_short(10);
57 FTList<FTSegment *> mid_short(10);
58 int n = _segments.length();
59 for (int i = 0; i^n; i++){
60 FTSegment * s = _segments[i];
61#ifndef OnlineMode
62 // s->printout();
63#endif
64 int stat = s->examine();
65 switch(stat){
66 case 0:
67 // simple
68 break;
69 case 1:
70 // inner short simple
71 inner_short.append(s);
72 n = _segments.remove(i);
73 break;
74 case 2:
75 // outer short simple
76 outer_short.append(s);
77 n = _segments.remove(i);
78 break;
79 case 3:
80 // to be divided
81 //if (!_superLayerId) _complecated_segments->append(s);
82 //else delete s;
83 _complecated_segments->append(s);
84 n = _segments.remove(i);
85 break;
86 default:
87 mid_short.append(s);
88 n = _segments.remove(i);
89 break;
90 }
91 }
92 connect_short_segments(inner_short,outer_short);
93 connect_singleHit(inner_short);
94 connect_singleHit(outer_short);
95 connect_singleHit(mid_short);
96 //if (!(_superLayerId&1)){
97 n = _segments.length();
98 for (int i = 0; i^n; i++){
99 _segments[i]->update();
100#ifndef OnlineMode
101 _segments[i]->printout();
102#endif
103 }
104}
105
106void
108{
109 // if(_superLayerId<5 || _superLayerId>9 ) return;
110 if( _superLayerId>9 ) return; // no need for the most outer SL
111
112 int n = _segments.length();
113 if(n>30) return;
114
115 float tdc[30][4]={0.},r[4]={0.},phi[4]={0.};
116 int wireId[4]={0};
117
118 IMdcGeomSvc* mdcGeomSvc;
119 StatusCode sc = Gaudi::svcLocator()->service("MdcGeomSvc", mdcGeomSvc);
120 // float size = mdcGeomSvc->Layer(0)->PCSiz();
121 float T1=0.5*0.1*(mdcGeomSvc->Layer(4*_superLayerId+0)->PCSiz())/0.004;
122 float T2=0.5*0.1*(mdcGeomSvc->Layer(4*_superLayerId+1)->PCSiz())/0.004;
123 float T3=0.5*0.1*(mdcGeomSvc->Layer(4*_superLayerId+2)->PCSiz())/0.004;
124 float T4=0.5*0.1*(mdcGeomSvc->Layer(4*_superLayerId+3)->PCSiz())/0.004;
125
126 for (int i = 0; i^n; i++){
127
128 float chi2=-99;
129 float t0c=0;
130
131 //float in_layerId=_segments[i]->innerBoundHits().first()->layer().layerId();
132 // float out_layerId=_segments[i]->outerBoundHits().first()->layer().layerId();
133 FTList<FTWire *>& hits=_segments[i]->wireHits();
134 // if(hits.length()>4) continue;
135
136 for(int j=0;j<hits.length();j++){
137 FTWire & h = *hits[j];
138
139 int layerId=h.layer().layerId();
140
141 /* if(wireId[layerId%4] == (h.localId()-1)||wireId[layerId%4]==(h.localId()+1)) {
142 tdc[i][layerId%4]=0;
143 break;
144 }*/
145
146 wireId[layerId%4]=h.localId();
147 phi[layerId%4]=h.phi();
148
149 // if(tdc[i][layerId%4] != (h.time())) continue;
150
151 tdc[i][layerId%4]=h.time();
152 r[layerId%4]=h.layer().r();
153 }
154 //// 4 cell
155
156 if(tdc[i][0]!=0 && tdc[i][1]!=0 && tdc[i][2]!=0 && tdc[i][3]!=0 && (wireId[0]==wireId[1]&& wireId[0]==wireId[2] && wireId[0]==wireId[3] || wireId[0]==wireId[1]-1 && wireId[0]==wireId[2]&& wireId[1]==wireId[3])){
157 // if(tdc[i][0]!=0 && tdc[i][1]!=0 && tdc[i][2]!=0 && tdc[i][3]!=0 &&( phi[0]<phi[1]&&phi[1]>phi[2]&&phi[2]<phi[3]||phi[0]>phi[1] && phi[1]<phi[2] && phi[2]>phi[3])){
158
159 float r0=r[0]-r[1]-(r[2]-r[1])/2;
160 float r1=-(r[2]-r[1])/2;
161 float r2=(r[2]-r[1])/2;
162 float r3=r[3]-r[2]+(r[2]-r[1])/2;
163
164 float t_i = tdc[i][0]+tdc[i][2];
165 float t_j = tdc[i][1]+tdc[i][3];
166 float l_j = T2+T4;
167 float r_i = r0+r2;
168 float r_j = r1+r3;
169 float r_2k= r0*r0+r1*r1+r2*r2+r3*r3;
170 float rt_i = r0*tdc[i][0]+r2*tdc[i][2];
171 float rt_j = r1*tdc[i][1]+r3*tdc[i][3];
172 float rl_j = r1*T2+r3*T4;
173
174 float deno= 4*r_2k-2*(r_i*r_i+r_j*r_j);
175
176 if (deno!=0){
177 float Pa=(4*(rt_i-rt_j+rl_j)-(t_i+t_j-l_j)*(r_i-r_j)-(t_i-t_j+l_j)*(r_i+r_j))/deno;
178 float Pb= 0.25*(t_i-t_j+l_j-(r_i+r_j) * Pa);
179 float Ang=fabs(90.-fabs(atan(Pa)*180./3.141593));
180
181 t0c = -0.25*((r_i-r_j)*Pa-t_i-t_j+l_j);
182
183 float pa1=(tdc[i][0]-tdc[i][2])/(r0-r2);
184 float pa2=(T4-T2-tdc[i][3]+tdc[i][1])/(r3-r1);
185
186 //// cout<<"t0c1: "<<-(tdc[i][0]-r0 *pa1 - 0.25*(t_i-t_j+l_j-(r_i+r_j)*pa1))<<endl;
187 //// cout<<"t0c2: "<<T2-tdc[i][1]-r1 *pa2 - 0.25*(t_i-t_j+l_j-(r_i+r_j)*pa2)<<endl;
188 //// cout<<"t0c3: "<<-(tdc[i][2]-r2 *pa1 - 0.25*(t_i-t_j+l_j-(r_i+r_j)*pa1))<<endl;
189 //// cout<<"t0c4: "<<T4-tdc[i][3]-r3 *pa2 - 0.25*(t_i-t_j+l_j-(r_i+r_j)*pa2)<<endl;
190
191 //// cout<<"t0c1,2,3,4: "<<-(tdc[i][0]-r0 * Pa-Pb)<<","<<(T2-tdc[i][1]-r1 * Pa-Pb)<<","<<-(tdc[i][2]-r2 * Pa-Pb)<<","<<(T4-tdc[i][3]-r3 * Pa-Pb)<<endl;
192
193 chi2=(tdc[i][0]-t0c-r0 * Pa-Pb)*(tdc[i][0]-t0c-r0 * Pa-Pb)+(T2-tdc[i][1]+t0c-r1 * Pa-Pb)*(T2-tdc[i][1]+t0c-r1 * Pa-Pb)+(tdc[i][2]-t0c-r2 * Pa-Pb)*(tdc[i][2]-t0c-r2 * Pa-Pb) + (T4-tdc[i][3]+t0c-r3 * Pa-Pb)*(T4-tdc[i][3]+t0c-r3 * Pa-Pb);
194
195 }
196 }
197 ////3 cell
198
199 /* else if(tdc[i][0]!=0 && tdc[i][1]!=0 && tdc[i][2]!=0 && (wireId[0]==wireId[1]-1||wireId[0]==wireId[1]) && wireId[2]==wireId[0]) {
200
201 //// else if(tdc[i][0]!=0 && tdc[i][1]!=0 && tdc[i][2]!=0 &&(phi[0]<phi[1] && phi[2]<phi[1] ||phi[0]>phi[1] && phi[2]>phi[1])){
202
203 float pa=-((r[0]+r[2]-2*r[1])*(tdc[i][0]+tdc[i][2])-2*(r[0]-r[1])*tdc[i][0]-2*(r[2]-r[1])*tdc[i][2])/((r[0]-r[2])*(r[0]-r[2]));
204 //float pa=(tdc[i][2]-tdc[i][0])/(r[2]-r[0]);
205 float pb=0.25*(tdc[i][0]-2*tdc[i][1]+tdc[i][2]+2.*T2-(r[0]+r[2]-2*r[1])*pa);
206 float Ang=fabs(90.-fabs(atan(pa)*180./3.141593));
207
208 t0c=-((r[0]+r[2]-2*r[1])*pa + 3.*pb - tdc[i][0]+tdc[i][1]-tdc[i][2]-T2);
209 //t0c=0.5*(tdc[i][1]+tdc[i][2]-(r[1]-r[2])*pa-T2);
210 chi2=(tdc[i][0]-t0c-pa*(r[0]-r[1])-pb)*(tdc[i][0]-t0c-pa*(r[0]-r[1])-pb)+(T2-tdc[i][1]+t0c-pb)*(T2-tdc[i][1]+t0c-pb)+(tdc[i][2]-t0c-pa*(r[2]-r[1])-pb)*(tdc[i][2]-t0c-pa*(r[2]-r[1])-pb);
211 }
212
213 else if(tdc[i][1]!=0 && tdc[i][2]!=0 && tdc[i][3]!=0 && (wireId[1]==wireId[2]||wireId[1]==wireId[2]+1) && wireId[3]==wireId[1]){
214 //// else if(tdc[i][1]!=0 && tdc[i][2]!=0 && tdc[i][3]!=0 &&(phi[1]<phi[2] && phi[3]<phi[2] ||phi[1]>phi[2] && phi[3]>phi[2])){
215
216 float pa=-((r[1]+r[3]-2*r[2])*(tdc[i][1]+tdc[i][3])-2*(r[1]-r[2])*tdc[i][1]-2*(r[3]-r[2])*tdc[i][3])/((r[1]-r[3])*(r[1]-r[3]));
217 float pb=0.25*(tdc[i][1]-2*tdc[i][2]+tdc[i][3]+2.*T3-(r[1]+r[3]-2*r[2])*pa);
218
219 t0c=-((r[1]+r[3]-2*r[2])*pa + 3.*pb - tdc[i][1]+tdc[i][2]-tdc[i][3]-T3);
220 chi2=(tdc[i][1]-t0c-pa*(r[1]-r[2])-pb)*(tdc[i][1]-t0c-pa*(r[1]-r[2])-pb)+(T3-tdc[i][2]+t0c-pb)*(T3-tdc[i][2]+t0c-pb)+(tdc[i][3]-t0c-pa*(r[3]-r[2])-pb)*(tdc[i][3]-t0c-pa*(r[3]-r[2])-pb);
221
222 }
223 */
224
225 // g_segmentchi2[_superLayerId][i]=chi2;
226 // g_t0c[_superLayerId][i]=t0c;
227
228 // if(_superLayerId>4 && chi2>0){
229 //cout<<"superLayer,t0c,chi2: "<<_superLayerId<<"["<<i<<"] : "<<t0c<<" , "<<chi2<<endl;
230 // }
231 // if(chi2<700 && _superLayerId>4 && chi2>-1){
232 if(chi2<(param->_chi2_segfit) && _superLayerId>4 && chi2>-1){
233 Estime[_superLayerId].append(t0c);
234 }
235 // if(chi2>0 && chi2<2000) continue;
236 // if (chi2<1000) continue;
237 continue;
238
239 _complecated_segments->append(_segments[i]);
240 for(int k=0;k<4;k++){
241 tdc[i][k]=0;
242 }
243 n = _segments.remove(i);
244 }
245}
246
247//cluster hits to segments
248void
249FTSuperLayer::clustering(void)
250{
251 // +---+---+
252 // | 4 | 5 |
253 // +-+-+---+-+-+ r
254 // Neighbor ID | 2 | * | 3 | ^
255 // +-+-+---+-+-+ |
256 // | 0 | 1 | +--> -phi
257 // +---+---+
258 //
259 if (!_wireHits.length()) return;
260 register FTWire ** hptr = _wireHits.firstPtr();
261 FTWire ** const last = _wireHits.lastPtr();
262 // discard continuous hits
263 do{(**hptr).chk_left_and_right();
264 }while((hptr++)!=last);
265
266 // clustering
267 hptr = _wireHits.firstPtr();
268 //FTList<FTWire *> * hits = new FTList<FTWire *>(10); zhangxm:because of 345 link
269 FTList<FTWire *> * hits = new FTList<FTWire *>(16);
270 do{ // loop over clusters
271 if ((**hptr).stateAND(FTWireHitAppendedOrInvalid)) continue;
272 int n = hits->append(*hptr);
273 (**hptr).stateOR(FTWireHitAppended);
274 for (int j = 0; j^n; j++){ // loop over hits in a cluster
275 const unsigned int checked = (*hits)[j]->state();
276 //const unsigned int * mptr = _neighborsMask;
277 FTWire** nptr = (*hits)[j]->neighborPtr();
278 // check 6 neighbors
279 for (unsigned Mask=FTWireNeighbor0; Mask^512; Mask<<=1, /*mptr++,*/ nptr++){
280 if ((**nptr).stateAND(FTWireHitAppendedOrInvalid)) continue;
281 n = hits->append(*nptr);
282 (**nptr).stateOR(FTWireHitAppended);
283 }
284 }
285
286 if (n^1){
287 _segments.append(new FTSegment(this, *hits));
288 //hits = new FTList<FTWire *>(10); because of 345 link
289 hits = new FTList<FTWire *>(16);
290 }else{
291 _singleHits.append(*hptr);
292 hits->clear();
293 }
294 }while((hptr++)!=last);
295 delete hits;
296}
297
298//to match two phi-near short segments and connect them
299void
300FTSuperLayer::connect_short_segments(FTList<FTSegment *> & inner_short,
301 FTList<FTSegment *> & outer_short)
302{
303 int n = inner_short.length();
304 int m = outer_short.length();
305 for (int i = 0; i^n; i++){
306 FTSegment * inner = inner_short[i];
307 const FTWire & in_outerBoundHit = * inner->outerBoundHits().first();
308 const FTLayer & in_outerBound = in_outerBoundHit.layer();
309 float in_outerPhi = inner->outgoingPhi();
310
311 int min_dphi_index = -1;
312 float min_dphi = 2*M_PI;
313 int dLayer_cache = -1;
314 float out_innerPhi, D_phi;
315 for (int j = 0; j^m; j++) {
316 const int out_innerLayerId
317 = (outer_short[j]->innerBoundHits().first())->layer().localLayerId();
318 int dLayer_cache_tmp = out_innerLayerId - in_outerBound.localLayerId();
319 if (dLayer_cache_tmp < 1) continue;
320 out_innerPhi = outer_short[j]->incomingPhi();
321 D_phi= fabs(in_outerPhi - out_innerPhi);
322 if (D_phi > M_PI) D_phi = 2*M_PI - D_phi;
323 if (D_phi < min_dphi) {
324 min_dphi = D_phi;
325 min_dphi_index = j;
326 dLayer_cache = dLayer_cache_tmp;
327 }
328 }
329 if (min_dphi_index < 0) continue;
330
331 switch (dLayer_cache) {
332 case 1:
333 if (min_dphi > _maxDphi[_superLayerId]*M_PI/180.) continue;
334 break;
335 case 2:
336 if (min_dphi > _maxDphi[_superLayerId]*M_PI/120.) continue;
337 break;
338 case 3:
339 if (min_dphi > _maxDphi[_superLayerId]*M_PI/80.) continue;
340 break;
341 default:
342 continue;
343 }
344
345 FTSegment * outer = outer_short[min_dphi_index];
346 n = inner_short.remove(i);
347 m = outer_short.remove(min_dphi_index);
348 inner->connect_outer(outer->wireHits(), outer->outerBoundHits());
349 _segments.append(inner);
350 delete outer;
351 }
352}
353
354
355void
356FTSuperLayer::connect_singleHit(FTList<FTSegment *> & short_sgmnts){
357 int minLength = 0;
358 if (superLayerId() == 2 || superLayerId() == 3 || superLayerId() == 4
359 || superLayerId() == 9 || superLayerId() == 10 )
360 minLength = 1;
361 int n = short_sgmnts.length();
362 int m = _singleHits.length();
363 for (int i = 0; i^n; i++){
364 FTSegment & s = * short_sgmnts[i];
365 const FTWire & outerBoundHit = * s.outerBoundHits().first();
366 const FTWire & innerBoundHit = * s.innerBoundHits().first();
367 const FTLayer & outerBound = outerBoundHit.layer();
368 const FTLayer & innerBound = innerBoundHit.layer();
369 float outerPhi = s.outgoingPhi();
370 float innerPhi = s.incomingPhi();
371
372 int min_dphi_index = -1;
373 float min_dphi = 2*M_PI;
374 int dLayer_cache = -1;
375 bool inner_flag_cache;
376 for (int j = 0; j^m; j++) {
377 const FTWire * h = _singleHits[j];
378 const int hLocalLayerId = h->layer().localLayerId();
379 int dLayer_cache_tmp = -1;
380 bool inner_flag_cache_tmp;
381 float D_phi;
382
383 if ((dLayer_cache_tmp=innerBound.localLayerId()-hLocalLayerId) > 0) {
384 D_phi = fabs(innerPhi - h->phi());
385 inner_flag_cache_tmp = true;
386 } else if ((dLayer_cache_tmp=hLocalLayerId-outerBound.localLayerId()) > 0) {
387 D_phi = fabs(outerPhi - h->phi());
388 inner_flag_cache_tmp = false;
389 } else {
390 continue;
391 }
392 if (D_phi > M_PI) D_phi = 2*M_PI - D_phi;
393
394 if (D_phi < min_dphi) {
395 min_dphi = D_phi;
396 min_dphi_index = j;
397 dLayer_cache = dLayer_cache_tmp;
398 inner_flag_cache = inner_flag_cache_tmp;
399 }
400 }
401
402 int not_append = 1;
403 if (min_dphi_index >= 0) {
404 switch (dLayer_cache) {
405 case 1:
406 if (min_dphi < _maxDphi[_superLayerId]*M_PI/180.) not_append = 0;
407 break;
408 case 2:
409 if (min_dphi < _maxDphi[_superLayerId]*M_PI/120.) not_append = 0;
410 break;
411 case 3:
412 if (min_dphi < _maxDphi[_superLayerId]*M_PI/80.) not_append = 0;
413 break;
414 default:
415 break;
416 }
417 }
418
419 if (!not_append) {
420 if (inner_flag_cache) {
421 s.connect_inner(_singleHits[min_dphi_index]);
422 } else {
423 s.connect_outer(_singleHits[min_dphi_index]);
424 }
425 _segments.append(&s);
426 m = _singleHits.remove(min_dphi_index);
427 } else {
428 if (outerBound.localLayerId()-innerBound.localLayerId()>minLength){
429 _segments.append(&s);
430 }else{
431 if (!_superLayerId) _complecated_segments->append(&s);
432 else delete &s;
433 }
434 }
435 }
436}
437
438
439void
441{
442 int n = _segments.length();
443 for (int i = 0; i^n; i++){
444 if (_segments[i]->track()) continue;
445 _complecated_segments->append(_segments[i]);
446 n = _segments.remove(i);
447 }
448}
449
450
451#ifdef FZISAN_DEBUG
452void
453FTSuperLayer::showBinary(unsigned src)
454{
455 unsigned mask = 512;
456 while(mask >>=1){
457 unsigned bit = (mask&src) ? 1 : 0;
458 }
459}
460#endif
const Int_t n
#define FTWireHitAppended
Definition FTWire.h:14
#define FTWireNeighbor0
Definition FTWire.h:16
#define FTWireHitInvalid
Definition FTWire.h:13
#define FTWireHitAppendedOrInvalid
Definition FTWire.h:15
XmlRpcServer s
#define M_PI
Definition TConstant.h:4
const float r(void) const
returns r form origin
Definition FTLayer.h:141
const int localLayerId(void) const
returns local-layer ID
Definition FTLayer.h:119
const int layerId(void) const
returns layer ID
Definition FTLayer.h:112
T & first(void) const
returns the first object in the list
Definition FTList.h:224
void deleteAll(void)
delete all object and clear(allocated memory remains same)
Definition FTList.h:266
void clear(void)
clear lists but the allocated memory remains same
Definition FTList.h:182
int remove(int &)
remove objects by index and returns decremented index and length
Definition FTList.h:140
T * lastPtr(void) const
returns the pointer of last object
Definition FTList.h:238
int length(void) const
returns the length of the list
Definition FTList.h:245
int append(const T &x)
append an object into the end of the list
Definition FTList.h:127
T * firstPtr(void) const
returns the pointer of first object
Definition FTList.h:231
float outgoingPhi(void) const
returns phi of outgoing position
Definition FTSegment.h:280
FTList< FTWire * > & wireHits(void) const
returns wire-hit list
Definition FTSegment.h:215
void connect_outer(const FTList< FTWire * > &, const FTList< FTWire * > &)
connect short segments
Definition FTSegment.h:170
FTList< FTWire * > & outerBoundHits(void) const
returns outerBoundHits
Definition FTSegment.h:229
static const float _maxDphi[11]
static MdcParameter * param
const int superLayerId(void) const
returns super-layer ID
void mkSegmentList(void)
create segment lists
void reduce_noise(FTList< float >(&Estime)[10])
calculate t0 and Chi2 from segment-fit in every superlayer
void reAppendSalvage(void)
append segments which are not used for tracks to the list for salvage
void clear(void)
clear object
float phi(void) const
returns phi
Definition FTWire.h:394
FTWire ** neighborPtr(void)
returns pointer of neighbor array
Definition FTWire.h:387
const FTLayer & layer(void) const
returns layer
Definition FTWire.h:359
const int localId(void) const
returns local ID
Definition FTWire.h:380
float time(void) const
rerurns TDC time(after t0 subtraction)
Definition FTWire.h:457
virtual const MdcGeoLayer *const Layer(unsigned id)=0
static MdcParameter * instance()
const float _chi2_segfit
Index first(Pair i)