00001
00002
00003
00004
00005
00006
00007
00008 #include "GaudiKernel/StatusCode.h"
00009 #include "GaudiKernel/Bootstrap.h"
00010
00011 #include "MdcGeomSvc/IMdcGeomSvc.h"
00012 #include "MdcGeomSvc/MdcGeoLayer.h"
00013
00014 #include "MdcFastTrkAlg/FTList.h"
00015 #include "MdcFastTrkAlg/FTLayer.h"
00016 #include "MdcFastTrkAlg/FTSuperLayer.h"
00017 #include "MdcFastTrkAlg/FTSegment.h"
00018 #include "MdcFastTrkAlg/FTWire.h"
00019 #include "MdcFastTrkAlg/MdcParameter.h"
00020 #include <iostream>
00021
00022 MdcParameter* FTSuperLayer::param = MdcParameter::instance();
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 const float FTSuperLayer::_maxDphi[11] =
00034 {11.9, 9.72, 6.26, 4.86, 3.81, 2.37, 2.08, 1.76, 1.53, 2.33, 1.88};
00035
00036 void
00037 FTSuperLayer::clear(void)
00038 {
00039 if (_wireHits.length()){
00040 register FTWire ** hptr = _wireHits.firstPtr();
00041 FTWire ** const last = _wireHits.lastPtr();
00042 do{(**hptr).state(FTWireHitInvalid);}while((hptr++)!=last);
00043 _wireHits.clear();
00044 }
00045 _singleHits.clear();
00046 if (_segments.length()) _segments.deleteAll();
00047 if (_complecated_segments) _complecated_segments->deleteAll();
00048 }
00049
00050
00051 void
00052 FTSuperLayer::mkSegmentList(void)
00053 {
00054 clustering();
00055 FTList<FTSegment *> inner_short(10);
00056 FTList<FTSegment *> outer_short(10);
00057 FTList<FTSegment *> mid_short(10);
00058 int n = _segments.length();
00059 for (int i = 0; i^n; i++){
00060 FTSegment * s = _segments[i];
00061 #ifndef OnlineMode
00062
00063 #endif
00064 int stat = s->examine();
00065 switch(stat){
00066 case 0:
00067
00068 break;
00069 case 1:
00070
00071 inner_short.append(s);
00072 n = _segments.remove(i);
00073 break;
00074 case 2:
00075
00076 outer_short.append(s);
00077 n = _segments.remove(i);
00078 break;
00079 case 3:
00080
00081
00082
00083 _complecated_segments->append(s);
00084 n = _segments.remove(i);
00085 break;
00086 default:
00087 mid_short.append(s);
00088 n = _segments.remove(i);
00089 break;
00090 }
00091 }
00092 connect_short_segments(inner_short,outer_short);
00093 connect_singleHit(inner_short);
00094 connect_singleHit(outer_short);
00095 connect_singleHit(mid_short);
00096
00097 n = _segments.length();
00098 for (int i = 0; i^n; i++){
00099 _segments[i]->update();
00100 #ifndef OnlineMode
00101 _segments[i]->printout();
00102 #endif
00103 }
00104 }
00105
00106 void
00107 FTSuperLayer::reduce_noise(FTList<float> (&Estime)[10])
00108 {
00109
00110 if( _superLayerId>9 ) return;
00111
00112 int n = _segments.length();
00113 if(n>30) return;
00114
00115 float tdc[30][4]={0.},r[4]={0.},phi[4]={0.};
00116 int wireId[4]={0};
00117
00118 IMdcGeomSvc* mdcGeomSvc;
00119 StatusCode sc = Gaudi::svcLocator()->service("MdcGeomSvc", mdcGeomSvc);
00120
00121 float T1=0.5*0.1*(mdcGeomSvc->Layer(4*_superLayerId+0)->PCSiz())/0.004;
00122 float T2=0.5*0.1*(mdcGeomSvc->Layer(4*_superLayerId+1)->PCSiz())/0.004;
00123 float T3=0.5*0.1*(mdcGeomSvc->Layer(4*_superLayerId+2)->PCSiz())/0.004;
00124 float T4=0.5*0.1*(mdcGeomSvc->Layer(4*_superLayerId+3)->PCSiz())/0.004;
00125
00126 for (int i = 0; i^n; i++){
00127
00128 float chi2=-99;
00129 float t0c=0;
00130
00131
00132
00133 FTList<FTWire *>& hits=_segments[i]->wireHits();
00134
00135
00136 for(int j=0;j<hits.length();j++){
00137 FTWire & h = *hits[j];
00138
00139 int layerId=h.layer().layerId();
00140
00141
00142
00143
00144
00145
00146 wireId[layerId%4]=h.localId();
00147 phi[layerId%4]=h.phi();
00148
00149
00150
00151 tdc[i][layerId%4]=h.time();
00152 r[layerId%4]=h.layer().r();
00153 }
00155
00156 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])){
00157
00158
00159 float r0=r[0]-r[1]-(r[2]-r[1])/2;
00160 float r1=-(r[2]-r[1])/2;
00161 float r2=(r[2]-r[1])/2;
00162 float r3=r[3]-r[2]+(r[2]-r[1])/2;
00163
00164 float t_i = tdc[i][0]+tdc[i][2];
00165 float t_j = tdc[i][1]+tdc[i][3];
00166 float l_j = T2+T4;
00167 float r_i = r0+r2;
00168 float r_j = r1+r3;
00169 float r_2k= r0*r0+r1*r1+r2*r2+r3*r3;
00170 float rt_i = r0*tdc[i][0]+r2*tdc[i][2];
00171 float rt_j = r1*tdc[i][1]+r3*tdc[i][3];
00172 float rl_j = r1*T2+r3*T4;
00173
00174 float deno= 4*r_2k-2*(r_i*r_i+r_j*r_j);
00175
00176 if (deno!=0){
00177 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;
00178 float Pb= 0.25*(t_i-t_j+l_j-(r_i+r_j) * Pa);
00179 float Ang=fabs(90.-fabs(atan(Pa)*180./3.141593));
00180
00181 t0c = -0.25*((r_i-r_j)*Pa-t_i-t_j+l_j);
00182
00183 float pa1=(tdc[i][0]-tdc[i][2])/(r0-r2);
00184 float pa2=(T4-T2-tdc[i][3]+tdc[i][1])/(r3-r1);
00185
00190
00192
00193 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);
00194
00195 }
00196 }
00198
00199
00200
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232 if(chi2<(param->_chi2_segfit) && _superLayerId>4 && chi2>-1){
00233 Estime[_superLayerId].append(t0c);
00234 }
00235
00236
00237 continue;
00238
00239 _complecated_segments->append(_segments[i]);
00240 for(int k=0;k<4;k++){
00241 tdc[i][k]=0;
00242 }
00243 n = _segments.remove(i);
00244 }
00245 }
00246
00247
00248 void
00249 FTSuperLayer::clustering(void)
00250 {
00251
00252
00253
00254
00255
00256
00257
00258
00259 if (!_wireHits.length()) return;
00260 register FTWire ** hptr = _wireHits.firstPtr();
00261 FTWire ** const last = _wireHits.lastPtr();
00262
00263 do{(**hptr).chk_left_and_right();
00264 }while((hptr++)!=last);
00265
00266
00267 hptr = _wireHits.firstPtr();
00268
00269 FTList<FTWire *> * hits = new FTList<FTWire *>(16);
00270 do{
00271 if ((**hptr).stateAND(FTWireHitAppendedOrInvalid)) continue;
00272 int n = hits->append(*hptr);
00273 (**hptr).stateOR(FTWireHitAppended);
00274 for (int j = 0; j^n; j++){
00275 const unsigned int checked = (*hits)[j]->state();
00276
00277 FTWire** nptr = (*hits)[j]->neighborPtr();
00278
00279 for (unsigned Mask=FTWireNeighbor0; Mask^512; Mask<<=1, nptr++){
00280 if ((**nptr).stateAND(FTWireHitAppendedOrInvalid)) continue;
00281 n = hits->append(*nptr);
00282 (**nptr).stateOR(FTWireHitAppended);
00283 }
00284 }
00285
00286 if (n^1){
00287 _segments.append(new FTSegment(this, *hits));
00288
00289 hits = new FTList<FTWire *>(16);
00290 }else{
00291 _singleHits.append(*hptr);
00292 hits->clear();
00293 }
00294 }while((hptr++)!=last);
00295 delete hits;
00296 }
00297
00298
00299 void
00300 FTSuperLayer::connect_short_segments(FTList<FTSegment *> & inner_short,
00301 FTList<FTSegment *> & outer_short)
00302 {
00303 int n = inner_short.length();
00304 int m = outer_short.length();
00305 for (int i = 0; i^n; i++){
00306 FTSegment * inner = inner_short[i];
00307 const FTWire & in_outerBoundHit = * inner->outerBoundHits().first();
00308 const FTLayer & in_outerBound = in_outerBoundHit.layer();
00309 float in_outerPhi = inner->outgoingPhi();
00310
00311 int min_dphi_index = -1;
00312 float min_dphi = 2*M_PI;
00313 int dLayer_cache = -1;
00314 float out_innerPhi, D_phi;
00315 for (int j = 0; j^m; j++) {
00316 const int out_innerLayerId
00317 = (outer_short[j]->innerBoundHits().first())->layer().localLayerId();
00318 int dLayer_cache_tmp = out_innerLayerId - in_outerBound.localLayerId();
00319 if (dLayer_cache_tmp < 1) continue;
00320 out_innerPhi = outer_short[j]->incomingPhi();
00321 D_phi= fabs(in_outerPhi - out_innerPhi);
00322 if (D_phi > M_PI) D_phi = 2*M_PI - D_phi;
00323 if (D_phi < min_dphi) {
00324 min_dphi = D_phi;
00325 min_dphi_index = j;
00326 dLayer_cache = dLayer_cache_tmp;
00327 }
00328 }
00329 if (min_dphi_index < 0) continue;
00330
00331 switch (dLayer_cache) {
00332 case 1:
00333 if (min_dphi > _maxDphi[_superLayerId]*M_PI/180.) continue;
00334 break;
00335 case 2:
00336 if (min_dphi > _maxDphi[_superLayerId]*M_PI/120.) continue;
00337 break;
00338 case 3:
00339 if (min_dphi > _maxDphi[_superLayerId]*M_PI/80.) continue;
00340 break;
00341 default:
00342 continue;
00343 }
00344
00345 FTSegment * outer = outer_short[min_dphi_index];
00346 n = inner_short.remove(i);
00347 m = outer_short.remove(min_dphi_index);
00348 inner->connect_outer(outer->wireHits(), outer->outerBoundHits());
00349 _segments.append(inner);
00350 delete outer;
00351 }
00352 }
00353
00354
00355 void
00356 FTSuperLayer::connect_singleHit(FTList<FTSegment *> & short_sgmnts){
00357 int minLength = 0;
00358 if (superLayerId() == 2 || superLayerId() == 3 || superLayerId() == 4
00359 || superLayerId() == 9 || superLayerId() == 10 )
00360 minLength = 1;
00361 int n = short_sgmnts.length();
00362 int m = _singleHits.length();
00363 for (int i = 0; i^n; i++){
00364 FTSegment & s = * short_sgmnts[i];
00365 const FTWire & outerBoundHit = * s.outerBoundHits().first();
00366 const FTWire & innerBoundHit = * s.innerBoundHits().first();
00367 const FTLayer & outerBound = outerBoundHit.layer();
00368 const FTLayer & innerBound = innerBoundHit.layer();
00369 float outerPhi = s.outgoingPhi();
00370 float innerPhi = s.incomingPhi();
00371
00372 int min_dphi_index = -1;
00373 float min_dphi = 2*M_PI;
00374 int dLayer_cache = -1;
00375 bool inner_flag_cache;
00376 for (int j = 0; j^m; j++) {
00377 const FTWire * h = _singleHits[j];
00378 const int hLocalLayerId = h->layer().localLayerId();
00379 int dLayer_cache_tmp = -1;
00380 bool inner_flag_cache_tmp;
00381 float D_phi;
00382
00383 if ((dLayer_cache_tmp=innerBound.localLayerId()-hLocalLayerId) > 0) {
00384 D_phi = fabs(innerPhi - h->phi());
00385 inner_flag_cache_tmp = true;
00386 } else if ((dLayer_cache_tmp=hLocalLayerId-outerBound.localLayerId()) > 0) {
00387 D_phi = fabs(outerPhi - h->phi());
00388 inner_flag_cache_tmp = false;
00389 } else {
00390 continue;
00391 }
00392 if (D_phi > M_PI) D_phi = 2*M_PI - D_phi;
00393
00394 if (D_phi < min_dphi) {
00395 min_dphi = D_phi;
00396 min_dphi_index = j;
00397 dLayer_cache = dLayer_cache_tmp;
00398 inner_flag_cache = inner_flag_cache_tmp;
00399 }
00400 }
00401
00402 int not_append = 1;
00403 if (min_dphi_index >= 0) {
00404 switch (dLayer_cache) {
00405 case 1:
00406 if (min_dphi < _maxDphi[_superLayerId]*M_PI/180.) not_append = 0;
00407 break;
00408 case 2:
00409 if (min_dphi < _maxDphi[_superLayerId]*M_PI/120.) not_append = 0;
00410 break;
00411 case 3:
00412 if (min_dphi < _maxDphi[_superLayerId]*M_PI/80.) not_append = 0;
00413 break;
00414 default:
00415 break;
00416 }
00417 }
00418
00419 if (!not_append) {
00420 if (inner_flag_cache) {
00421 s.connect_inner(_singleHits[min_dphi_index]);
00422 } else {
00423 s.connect_outer(_singleHits[min_dphi_index]);
00424 }
00425 _segments.append(&s);
00426 m = _singleHits.remove(min_dphi_index);
00427 } else {
00428 if (outerBound.localLayerId()-innerBound.localLayerId()>minLength){
00429 _segments.append(&s);
00430 }else{
00431 if (!_superLayerId) _complecated_segments->append(&s);
00432 else delete &s;
00433 }
00434 }
00435 }
00436 }
00437
00438
00439 void
00440 FTSuperLayer::reAppendSalvage(void)
00441 {
00442 int n = _segments.length();
00443 for (int i = 0; i^n; i++){
00444 if (_segments[i]->track()) continue;
00445 _complecated_segments->append(_segments[i]);
00446 n = _segments.remove(i);
00447 }
00448 }
00449
00450
00451 #ifdef FZISAN_DEBUG
00452 void
00453 FTSuperLayer::showBinary(unsigned src)
00454 {
00455 unsigned mask = 512;
00456 while(mask >>=1){
00457 unsigned bit = (mask&src) ? 1 : 0;
00458 }
00459 }
00460 #endif