00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include "MdcGeom/BesAngle.h"
00018 #include "MdcData/MdcHit.h"
00019 #include "CLHEP/Alist/AList.h"
00020 #include "MdcGeom/MdcDetector.h"
00021 #include "MdcGeom/MdcSuperLayer.h"
00022 #include "MdcTrkRecon/MdcSeg.h"
00023 #include "MdcTrkRecon/MdcSegGrouperSt.h"
00024 #include "MdcTrkRecon/MdcSegList.h"
00025 #include "MdcTrkRecon/MdcSegInfoSterO.h"
00026 #include "MdcTrkRecon/MdcSegWorks.h"
00027 #include "MdcTrkRecon/MdcSegUsage.h"
00028 #include "MdcTrkRecon/MdcTrack.h"
00029 #include "MdcTrkRecon/GmsList.h"
00030 #include "TrkBase/TrkExchangePar.h"
00031 #include "TrkBase/TrkRecoTrk.h"
00032 #include "TrkBase/TrkFit.h"
00033 #include "TrkFitter/TrkHelixMaker.h"
00034
00035
00036 #include "AIDA/IHistogram1D.h"
00037 extern AIDA::IHistogram1D* g_z0Cut;
00038 extern AIDA::IHistogram1D* g_ctCut;
00039 extern AIDA::IHistogram1D* g_delCt;
00040 extern AIDA::IHistogram1D* g_delZ0;
00041
00042
00043
00044
00045 MdcSegGrouperSt::MdcSegGrouperSt(const MdcDetector *gm, int debug) :
00046 MdcSegGrouper(gm, gm->nStereoSuper(-1) + gm->nStereoSuper(1), debug) {
00047
00048
00049 lTestGroup = true;
00050 lTestSingle = false;
00051
00052 isValid = new bool * [nPly()];
00053 for (int j = 0; j < nPly(); j++) {
00054 isValid[j] = 0;
00055 }
00056 }
00057
00058
00059 void MdcSegGrouperSt::resetList() {
00060
00061
00062 int nsuper = _gm->nSuper();
00063 for (int i = 0; i < nsuper; i++) {
00064 segList[i].removeAll();
00065 }
00066 }
00067
00068
00069 void MdcSegGrouperSt::fillWithSegs(const MdcSegList *inSegs,
00070 const MdcTrack *track) {
00071
00072
00073
00074
00075 int nsuper = _gm->nSuper();
00076 resetList();
00077 if(3==_debug) std::cout<<std::endl<< "=====MdcSegGrouperSt::fillWithSegs====="<<std::endl;
00078
00079 const TrkFit* tkFit = track->track().fitResult();
00080 if (tkFit == 0) return;
00081 TrkExchangePar par = tkFit->helix(0.0);
00082 double kap = 0.5 * par.omega();
00083 double phi0 = par.phi0();
00084 double d0 = par.d0();
00085 MdcSegWorks segStuff;
00086
00087 bool lStraight = false;
00088 if (fabs(kap) < 1.e-9) { lStraight = true; }
00089
00090 double rCurl = 99999999.;
00091 if (!lStraight) {
00092 rCurl = fabs(d0 + 1./kap);
00093 }
00094
00095
00096 MdcSegInfoSterO *info = new MdcSegInfoSterO;
00097
00098
00099 for (int isuper = 0; isuper < nsuper; isuper++) {
00100 const GmsList *inList = inSegs->oneList(isuper);
00101
00102 if (inList->count() == 0) continue;
00103
00104 MdcSeg *inSeg = (MdcSeg *) inList->first();
00105
00106 if (inSeg->superlayer()->whichView() == 0) continue;
00107
00108
00109 double radius = inSeg->superlayer()->rEnd();
00110 if (radius >= rCurl) break;
00111 double phiArg = kap * radius;
00112 if (lStraight) phiArg = d0 / radius;
00113 if (phiArg < -1.0) phiArg = -1.0;
00114 else if (phiArg > 1.0) phiArg = 1.0;
00115 segStuff.phiArg = phiArg;
00116 segStuff.phiAx.setRad(phi0 + asin(phiArg));
00117 BesAngle delPhi( fabs(inSeg->superlayer()->delPhi()) );
00118
00119
00120 BesAngle maxPhiA = segStuff.phiAx + delPhi + 0.1;
00121 BesAngle minPhiA = segStuff.phiAx - delPhi - 0.1;
00122
00123 double maxPhi = maxPhiA;
00124 double minPhi = minPhiA;
00125 bool phitest = (maxPhi > minPhi);
00126
00127
00128
00129
00130 segStuff.wirLen2inv = 1./ (inSeg->superlayer()->zEnd() *
00131 inSeg->superlayer()->zEnd() );
00132
00133 if(_debug==3) std::cout<<std::endl<<"Pick segment by phi from " <<minPhi<< " to "<<maxPhi<<" phiAx="<<segStuff.phiAx<<" delPhi="<<delPhi<<std::endl;
00134 for ( ; inSeg != 0; inSeg = (MdcSeg *) inSeg->next()) {
00135
00136 BesAngle phiSeg(inSeg->phi());
00137
00138 if(_debug==3){ inSeg->plotSeg(); }
00139 if (phitest) {
00140 if (phiSeg < minPhi){
00141 if (_debug ==3){std::cout << " CUT by phi "<<phiSeg
00142 << "<min "<<minPhi << std::endl;}
00143 continue;
00144 }
00145 else if(phiSeg > maxPhi) {
00146 if (_debug ==3){std::cout << " CUT by phi "<<phiSeg
00147 <<">max "<<maxPhi<< std::endl;}
00148 break;
00149 }
00150 } else {
00151 if (phiSeg > maxPhi && phiSeg < minPhi) {
00152 if (_debug ==3){std::cout << "!phitest "<<phiSeg
00153 <<" max "<<maxPhi<< " min "<<minPhi << std::endl;}
00154 continue;
00155 }
00156 }
00157
00158 if(_debug == 3) std::cout<<" **KEEP seg phi="<<phiSeg<< std::endl;
00159
00160
00161 int isBad = info->calcStereo(inSeg, track->track(), segStuff);
00162
00163 if (isBad) {
00164 if (_debug ==3){std::cout << " CUT by calcStereo isBad!"<<std::endl;}
00165 continue;
00166 }
00167
00168 if (g_z0Cut) {g_z0Cut->fill(info->z0());}
00169 if (g_ctCut) {g_ctCut->fill(info->ct());}
00170 if (fabs(info->ct()) > inSeg->segPar()->ctcut) {
00171 if (_debug ==3){std::cout << " CUT by ctcut! "
00172 <<fabs(info->ct())<<" > cut:"<< inSeg->segPar()->ctcut<<" "<<phiSeg<<std::endl; }
00173 continue;
00174 }
00175 if (fabs(info->z0()) > inSeg->segPar()->z0cut){
00176 if (_debug ==3){std::cout << " CUT by z0cut! "
00177 <<fabs(info->z0())<<" >cut "<< inSeg->segPar()->z0cut<<" "<<phiSeg<<std::endl; }
00178 continue;
00179 }
00180 inSeg->setInfo(info);
00181 info = new MdcSegInfoSterO;
00182 segList[isuper].append(inSeg);
00183
00184
00185 }
00186 }
00187 delete info;
00188 }
00189
00190
00191 void MdcSegGrouperSt::plotStereo( ) const {
00192
00193 int nsuper = _gm->nSuper();
00194 int isuper;
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207 for (isuper = 0; isuper < nsuper; isuper++) {
00208 if (segList[isuper].length() == 0) continue;
00209
00210 MdcSeg *inSeg = (MdcSeg *) segList[isuper].first();
00211
00212
00213 for (int j = 0 ; j < (int) segList[isuper].length(); j++) {
00214 segList[isuper][j]->plotSeg();
00215 }
00216 }
00217
00218 }
00219
00220 int MdcSegGrouperSt::incompWithSeg(const MdcSeg *refSeg,
00221 const MdcSeg *testSeg) {
00222
00223
00224 return 0;
00225 }
00226
00227
00228 int MdcSegGrouperSt::incompWithGroup(MdcSeg **segGroup, const MdcSeg *testSeg,
00229 int iply) {
00230
00231
00232
00233
00234
00235
00236 int i;
00237
00238 for (i = iply - 1; i > 0; i--) {
00239 if (!leaveGap[i]) break;
00240 }
00241 const MdcSeg *first = (*combList[i])[currentSeg[i]];
00242
00243
00244 MdcSegInfoSterO* firstInfo = (MdcSegInfoSterO *) first->info();
00245 MdcSegInfoSterO* newInfo = (MdcSegInfoSterO *) testSeg->info();
00246
00247 if (g_delCt) {g_delCt->fill( firstInfo->ct() - newInfo->ct());}
00248 if (fabs( firstInfo->ct() - newInfo->ct() ) >
00249 testSeg->segPar()->delCtCut) {
00250 if(_debug==3){
00251 cout << "---MdcSegGrouperSt Ct CUT!" << endl;
00252 std::cout<<"first:"; first->plotSeg(); std::cout<<std::endl;
00253 std::cout<<"test:"; testSeg->plotSeg();std::cout<<std::endl;
00254 std::cout << "first.ct "<< firstInfo->ct()
00255 <<" test.ct "<<newInfo->ct()
00256 <<" delta ct "<<firstInfo->ct() - newInfo->ct()
00257 <<" Cut "<<testSeg->segPar()->delCtCut << std::endl;
00258
00259 std::cout<<"--- "<< std::endl;
00260 }
00261 return -1;
00262 }else{
00263 }
00264
00265 double arcFirst = firstInfo->arc();
00266 double arcNew = newInfo->arc();
00267 double zFirst = firstInfo->z0() + firstInfo->ct() * arcFirst;
00268 double zNew = newInfo->z0() + newInfo->ct() * arcNew;
00269
00270 double zProj = zFirst / arcFirst * arcNew;
00271 if (g_delZ0) {g_delZ0->fill( zProj - zNew);}
00272
00273 if (fabs(zProj - zNew) > testSeg->segPar()->delZ0Cut) {
00274 if(3==_debug){
00275 cout << "MdcSegGrouperSt delZ0Cut not incompWithGroup CUT!" << endl;
00276 testSeg->plotSeg();
00277 std::cout<<" zProj "<< zProj << " zNew "<< zNew<< " CUT! "
00278 << testSeg->segPar()->delZ0Cut << std::endl;
00279 }
00280 return -1;
00281 }
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314 return 0;
00315 }
00316
00317
00318 void MdcSegGrouperSt::resetComb(const MdcSeg *seed) {
00319
00320
00321
00322 if (isValid != 0) {
00323 int i;
00324 for (i = 0; i < nDeep; i++) {
00325 delete [] isValid[i];
00326 isValid[i] = 0;
00327 }
00328 }
00329
00330
00331 int islay = 0;
00332 int iply = 0;
00333 nPlyFilled = 0;
00334 nNull = 0;
00335
00336
00337
00338 for (const MdcSuperLayer *thisSlay = _gm->firstSlay(); thisSlay != 0;
00339 thisSlay = thisSlay->next()) {
00340
00341 islay++;
00342 if (thisSlay->whichView() == 0) continue;
00343 firstGood[iply] = 0;
00344 firstBad[iply] = segList[islay-1].length();
00345 if (firstGood[iply] > firstBad[iply]) firstGood[iply] = firstBad[iply];
00346 if (firstGood[iply] == firstBad[iply]) {
00347
00348 continue;
00349 }
00350
00351 combList[iply] = &segList[islay-1];
00352 leaveGap[iply] = false;
00353 iply++;
00354 }
00355
00356 nPlyFilled = iply;
00357 resetSegCounters();
00358 maxNull = nPlyFilled - 2;
00359 }
00360
00361
00362 MdcTrack* MdcSegGrouperSt::storePar(MdcTrack* trk, double parms[2],
00363 double chisq,
00364 TrkContext&, double) {
00365
00366 assert (trk != 0);
00367 TrkHelixMaker maker;
00368 if(3==_debug) std::cout<< "-----storePar z0 "<<parms[0]<<" ct "<<parms[1]<<std::endl;
00369 maker.addZValues(trk->track(), parms[0], parms[1], chisq);
00370 return trk;
00371 }
00372