00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include "MdcGeom/Constants.h"
00018 #include "MdcGeom/BesAngle.h"
00019 #include "MdcData/MdcHit.h"
00020 #include "MdcData/MdcHitUse.h"
00021 #include "MdcGeom/MdcDetector.h"
00022 #include "MdcGeom/MdcSuperLayer.h"
00023 #include "MdcTrkRecon/MdcSegList.h"
00024 #include "MdcTrkRecon/mdcWrapWire.h"
00025 #include "MdcTrkRecon/MdcSeg.h"
00026 #include "MdcTrkRecon/GmsList.h"
00027
00028 #include "AIDA/IHistogram1D.h"
00029 extern AIDA::IHistogram1D* g_phiDiff;
00030 extern AIDA::IHistogram1D* g_slopeDiff;
00031
00032 const double _twopi = Constants::twoPi;
00033
00034
00035 MdcSegList::MdcSegList(int nSlayer, const MdcSegParams segP) {
00036
00037 segParam = segP;
00038 MdcSeg::setParam(&segParam);
00039 zeroSeed();
00040 segList = new GmsList[nSlayer];
00041 _nsupers = nSlayer;
00042 }
00043
00044
00045 MdcSegList::~MdcSegList() {
00046
00047 delete [] segList;
00048 }
00049
00050
00051 void MdcSegList::zeroSeed() {
00052
00053 seedSeg[0] = seedSeg[1] = seedSeg[2] = 0;
00054 seedSlay[0] = seedSlay[1] = seedSlay[2] = 0;
00055 }
00056
00057 void MdcSegList::resetSeed(const MdcDetector *gm) {
00058
00059
00060 for (int iview = -1; iview <= 1; iview++) {
00061 int viewIndex = iview + 1;
00062
00063
00064
00065
00066
00067
00068
00069 seedSlay[viewIndex] = gm->lastSlayInView(iview);
00070 if (seedSlay[viewIndex] != 0) seedSeg[viewIndex] =
00071 (MdcSeg *) oneList(seedSlay[viewIndex]->index())->first();
00072 }
00073 }
00074
00075 void MdcSegList::plot() const {
00076
00077
00078 MdcSeg *aSeg;
00079 for (int i = 0; i < _nsupers; i++) {
00080 std::cout << " ---------------MdcSeg of Slayer "<<i<<"-------------" << std::endl;
00081
00082 aSeg = (MdcSeg *) segList[i].first();
00083 while (aSeg != 0) {
00084 aSeg->plotSegAll();
00085 std::cout << endl;
00086 aSeg = (MdcSeg *) aSeg->next();
00087 }
00088 }
00089 }
00090
00091
00092 void MdcSegList::destroySegs() {
00093
00094
00095 zeroSeed();
00096 MdcSeg *aSeg, *bSeg;
00097 for (int i = 0; i < _nsupers; i++) {
00098 aSeg = (MdcSeg *) segList[i].first();
00099 while (aSeg != 0) {
00100 bSeg = (MdcSeg *) aSeg->next();
00101 segList[i].remove( aSeg );
00102 delete aSeg ;
00103 aSeg = bSeg;
00104 }
00105 }
00106 return;
00107 }
00108
00109
00110 void MdcSegList::sortByPhi() {
00111
00112
00113
00114 for (int isuper = 0; isuper < _nsupers; isuper++) {
00115
00116 if (segList[isuper].first() == 0) continue;
00117 MdcSeg *aseg = (MdcSeg *) segList[isuper].first()->next();
00118 for (int j = 1; j < (int) segList[isuper].count(); j++) {
00119 double phiSeg = aseg->phi();
00120 MdcSeg *nextseg = (MdcSeg *) aseg->next();
00121
00122 MdcSeg *bseg = (MdcSeg *) aseg->prev();
00123 if (phiSeg < bseg->phi() ) {
00124 while (bseg != 0 && phiSeg < bseg->phi()) {
00125 bseg = (MdcSeg *) bseg->prev();
00126 }
00127 segList[isuper].moveAfter(aseg, bseg);
00128 }
00129
00130 aseg = nextseg;
00131 }
00132
00133 }
00134 return;
00135 }
00136
00137
00138
00139 MdcSeg * MdcSegList::getSeed( int iview, int goodOnly) {
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149 int viewIndex = iview + 1;
00150 if (seedSlay[viewIndex] == 0) return 0;
00151 while (1) {
00152
00153
00154
00155
00156 if (seedSeg[viewIndex] == 0) {
00157
00158 seedSlay[viewIndex] = seedSlay[viewIndex]->prevInView();
00159 if (seedSlay[viewIndex] == 0) return 0;
00160 seedSeg[viewIndex] = (MdcSeg *)
00161 oneList(seedSlay[viewIndex]->index())->first();
00162 }
00163 else if (seedSeg[viewIndex]->segUsed()) {
00164
00165 seedSeg[viewIndex] = (MdcSeg *) seedSeg[viewIndex]->next();
00166 }
00167 else if (seedSeg[viewIndex]->segAmbig() && goodOnly) {
00168
00169 seedSeg[viewIndex] = (MdcSeg *) seedSeg[viewIndex]->next();
00170
00171
00172 }
00173 else if ( !seedSeg[viewIndex]->segFull() && goodOnly) {
00174
00175 seedSeg[viewIndex] = (MdcSeg *) seedSeg[viewIndex]->next();
00176
00177 } else {
00178
00179 int nused = 0;
00180 int navail = seedSeg[viewIndex]->nHit();
00181 for (int ihit = 0; ihit < seedSeg[viewIndex]->nHit(); ihit++) {
00182 MdcHitUse *ahit = seedSeg[viewIndex]->hit(ihit);
00183 if ( ahit->mdcHit()->usedHit() ) nused++;
00184 }
00185
00186 if (navail - nused < 2 && navail > 0) {
00187 seedSeg[viewIndex] = (MdcSeg *) seedSeg[viewIndex]->next();
00188 } else {
00189
00190 break;
00191 }
00192 }
00193
00194
00195
00196
00197
00198
00199
00200 }
00201
00202
00203 MdcSeg *theSeed = seedSeg[viewIndex];
00204 seedSeg[viewIndex] = (MdcSeg *) seedSeg[viewIndex]->next();
00205
00206 return theSeed;
00207 }
00208
00209 void MdcSegList::tagAmbig() {
00210
00211
00212
00213
00214 for (int isuper = 0; isuper < _nsupers; isuper++) {
00215 const GmsList *thisSegList = &segList[isuper];
00216 if (thisSegList->count() < 2) return;
00217
00218
00219 double width, diff;
00220 double thisphi, nextphi;
00221 double nisocell = 2;
00222
00223
00224
00225
00226 const MdcLayer *aLayer =
00227 ((MdcSeg *)thisSegList->first())->superlayer()->firstLayer();
00228 width = _twopi / aLayer->nWires();
00229 diff = nisocell * width;
00230
00231
00232 for ( MdcSeg *startSeg = (MdcSeg *) thisSegList->first();
00233 startSeg != 0; startSeg = (MdcSeg *) startSeg->next() ) {
00234 MdcSeg *nextSeg = (MdcSeg *) startSeg->next();
00235 if (nextSeg == 0) nextSeg = (MdcSeg *) thisSegList->first();
00236
00237 thisphi = startSeg->phi();
00238 nextphi = nextSeg->phi();
00239 if (nextphi < thisphi) nextphi += _twopi;
00240
00241 if (nextphi - thisphi < diff) {
00242 startSeg->setAmbig();
00243 nextSeg->setAmbig();
00244 }
00245
00246 }
00247 }
00248
00249 return;
00250 }
00251
00252
00253
00254
00255
00256
00257 void MdcSegList::append(MdcSeg *aSeg) {
00258
00259 segList[aSeg->superlayer()->index()].append(aSeg);
00260 }
00261
00262 void MdcSegList::massageSegs() {
00263
00264 sortByPhi();
00265
00266 bool drop = (segPar()->dropDups && count() > 200);
00267 deleteDups(drop);
00268 tagAmbig();
00269 }
00270
00271 void MdcSegList::deleteDups(bool ldrop) {
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283 bool printit = (5 == segPar()->lPrint);
00284
00285
00286 if(printit) std::cout<< " =======MdcSegList::deleteDups()===== " << std::endl;
00287
00288
00289 for (int isuper = 0; isuper < _nsupers; isuper++) {
00290 if(isuper==10) continue;
00291 GmsList *thisSegList = &segList[isuper];
00292 if (thisSegList->count() < 2) continue;
00293 if(5 == segPar()->lPrint || 13 == segPar()->lPrint){
00294 std::cout<<endl<< " -----slayer " << isuper<<"-----"<< std::endl;
00295 }
00296
00297 double thisphi, nextphi;
00298 double slopediff = segPar()->slopeDiffDrop;
00299
00300
00301 const MdcLayer *aLayer =
00302 ((MdcSeg *)thisSegList->first())->superlayer()->firstLayer();
00303 double width = _twopi / aLayer->nWires();
00304 double phidiff = segPar()->phiDiffDropMult * width;
00305
00306 int lWrapped = 0;
00307
00308 MdcSeg *startSeg = (MdcSeg *) thisSegList->first();
00309 while ( startSeg != 0 ) {
00310 if (thisSegList->count() < 2) break;
00311 if (lWrapped) break;
00312 thisphi = startSeg->phi();
00313 MdcSeg *nextSeg = startSeg;
00314 if(5 == segPar()->lPrint || 13 == segPar()->lPrint){
00315 startSeg->plotSeg();std::cout<<std::endl;
00316 }
00317
00318 while (1) {
00319 nextSeg = (MdcSeg *) nextSeg->next();
00320 if (nextSeg == 0) {
00321 nextSeg = (MdcSeg *) thisSegList->first();
00322 if (lWrapped == 1) break;
00323 lWrapped = 1;
00324 }
00325 if( nextSeg == startSeg ) break;
00326 nextphi = nextSeg->phi();
00327
00328 if (nextphi < thisphi) nextphi += _twopi;
00329
00330 if(printit){
00331 std::cout<<std::endl<<"----start "; startSeg->plotSeg(); std::cout<<std::endl;
00332 std::cout<<" next "; nextSeg->plotSeg(); std::cout<<std::endl;
00333 std::cout<<" phi diff "<< nextphi -thisphi << " cut phidiff "<<phidiff;
00334 std::cout<<" slope diff "<< nextSeg->slope() - startSeg->slope()
00335 << " cut slopdiff "<<slopediff<<std::endl;
00336 }
00337
00338 if (g_phiDiff) {g_phiDiff->fill((nextphi - thisphi)/width);}
00339 if (g_slopeDiff) {g_slopeDiff->fill(nextSeg->slope() - startSeg->slope());}
00340
00341 if (nextphi - thisphi > phidiff) {
00342 if(printit){ std::cout<< " --not same phi" <<std::endl; }
00343
00344 startSeg = (MdcSeg *) startSeg->next();
00345 break;
00346 } else if (fabs( nextSeg->slope() - startSeg->slope() ) > slopediff) {
00347 if(printit){ std::cout<< " --attached, but not identical -- move on to next nextSeg "<< std::endl; }
00348
00349 continue;
00350 } else {
00351
00352 int firstBetter = 0;
00353 assert (startSeg->nHit() > 2);
00354 double chiFirst = startSeg->chisq() / (startSeg->nHit() - 2);
00355 assert (nextSeg->nHit() > 2);
00356 double chiNext = nextSeg->chisq() / (nextSeg->nHit() - 2);
00357 double theDiff;
00358 int cdiff = (int) nextSeg->nHit() - (int) startSeg->nHit();
00359 theDiff = (chiFirst - chiNext) + 2. * (cdiff);
00360 if (theDiff < 0.) firstBetter = 1;
00361 if(printit){
00362 std::cout<< " --identical; choose better "<< std::endl;
00363 std::cout<< " chi start "<<chiFirst<<" chi next "<< chiNext<<" cdiff "<<cdiff<< std::endl;
00364 }
00365
00366 if (firstBetter) {
00367 if(printit){ std::cout<< " startSeg better,"; }
00368 MdcSeg *tempSeg = nextSeg;
00369 nextSeg = (MdcSeg *) nextSeg->next();
00370 if (ldrop) {
00371 if(printit){ std::cout<< " delete nextSeg"<<std::endl; }
00372 thisSegList->remove(tempSeg);
00373 delete tempSeg;
00374 } else {
00375 if(printit){ std::cout<< " nextSeg->setUsed(true) "<<std::endl; }
00376 tempSeg->setUsed();
00377 }
00378 if (nextSeg == 0) {
00379 nextSeg = (MdcSeg *) thisSegList->first();
00380 if (lWrapped == 1) break;
00381 lWrapped = 1;
00382 }
00383 } else {
00384 if(printit){ std::cout<< " nextSeg better,"; }
00385 MdcSeg *tempSeg = startSeg;
00386 startSeg = (MdcSeg *) startSeg->next();
00387 if (ldrop) {
00388 if(printit){ std::cout<< " delete startSeg"<<std::endl; }
00389 thisSegList->remove(tempSeg);
00390 delete tempSeg;
00391 } else {
00392 if(printit){ std::cout<< " startSeg->setUsed(true) "<<std::endl; }
00393 tempSeg->setUsed();
00394 }
00395 break;
00396 }
00397
00398 }
00399 }
00400
00401 }
00402 }
00403 }
00404
00405 void
00406 MdcSegList::setPlot(int lPlt) {
00407 MdcSeg::segPar()->lPlot = lPlt;
00408 }
00409
00410
00411 const GmsList*
00412 MdcSegList::oneList(int slayIndex) const {
00413
00414 return &segList[slayIndex];
00415 }
00416
00417
00418 int
00419 MdcSegList::count() const {
00420
00421 int cnt = 0;
00422 for (int i = 0; i < _nsupers; i++) cnt += segList[i].count();
00423 return cnt;
00424 }