00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include "MdcGeom/BesAngle.h"
00018 #include "CLHEP/Alist/AList.h"
00019 #include "MdcData/MdcHit.h"
00020 #include "MdcGeom/MdcDetector.h"
00021 #include "MdcGeom/MdcSuperLayer.h"
00022 #include "MdcTrkRecon/MdcSegGrouperAx.h"
00023 #include "MdcTrkRecon/MdcSegList.h"
00024 #include "MdcTrkRecon/MdcSegInfoAxialO.h"
00025 #include "MdcTrkRecon/MdcSegUsage.h"
00026 #include "MdcTrkRecon/MdcTrack.h"
00027 #include "MdcTrkRecon/mdcWrapAng.h"
00028 #include "MdcTrkRecon/MdcSeg.h"
00029 #include "MdcTrkRecon/GmsList.h"
00030 #include "TrkBase/TrkExchangePar.h"
00031 extern double MdcTrkReconCut_combAxPhi0;
00032 extern double MdcTrkReconCut_combAxCurv;
00033 extern double MdcTrkReconCut_combAxPhi0Cut;
00034 extern double MdcTrkReconCut_combAxCurvCut;
00035
00036
00037 #include "GaudiKernel/NTuple.h"
00038 extern NTuple::Tuple* g_tupleCombAx;
00039 extern NTuple::Item<double> g_combAxdPhi0;
00040 extern NTuple::Item<double> g_combAxdCurv;
00041 extern NTuple::Item<double> g_combAxSigPhi0;
00042 extern NTuple::Item<double> g_combAxSigCurv;
00043 extern NTuple::Item<double> g_combAxSlSeed;
00044 extern NTuple::Item<double> g_combAxSlTest;
00045 extern NTuple::Item<double> g_combAxQualitySeed;
00046 extern NTuple::Item<double> g_combAxQualityTest;
00047 extern NTuple::Item<double> g_combAxPatSeed;
00048 extern NTuple::Item<double> g_combAxPatTest;
00049 extern NTuple::Item<double> g_combAxNhitSeed;
00050 extern NTuple::Item<double> g_combAxNhitTest;
00051 extern NTuple::Item<double> g_combAxMc;
00052 extern NTuple::Item<double> g_combAxMcPt;
00053 extern NTuple::Item<double> g_combAxMcTheta;
00054 extern NTuple::Item<double> g_combAxMcPhi;
00055 extern NTuple::Item<double> g_combAxMcAmbigSeed;
00056 extern NTuple::Item<double> g_combAxMcAmbigTest;
00057
00058
00059
00060 MdcSegGrouperAx::MdcSegGrouperAx(const MdcDetector *gm, int debug) :
00061 MdcSegGrouper(gm, gm->nAxialSuper() - 1, debug) {
00062
00063 lTestGroup = false;
00064 lTestSingle = true;
00065
00066 isValid = new bool * [nPly()];
00067 for (int j = 0; j < nPly(); j++) {
00068 isValid[j] = 0;
00069 }
00070 }
00071
00072
00073 void MdcSegGrouperAx::fillWithSegs( const MdcSegList *inSegs) {
00074
00075 if(3==_debug) std::cout<<std::endl<< "=====MdcSegGrouperAx::fillWithSegs====="<<std::endl;
00076
00077
00078 for (int isuper = 0; isuper < _gm->nSuper(); isuper++) {
00079 const GmsList *inList = inSegs->oneList(isuper);
00080 if (inList->count() == 0) continue;
00081 MdcSeg *inSeg = (MdcSeg *) inList->first();
00082
00083 if (inSeg->superlayer()->whichView() != 0) continue;
00084
00085 while (inSeg != 0) {
00086
00087 MdcSegInfoAxialO *info = new MdcSegInfoAxialO;
00088 inSeg->setInfo(info);
00089 info->calcFromOrigin(inSeg);
00090
00091
00092
00093 int isInserted = 0;
00094 for (int iseg = 0; iseg < (int) segList[isuper].length(); iseg++) {
00095 MdcSeg *aSeg = segList[isuper][iseg];
00096 if ( ((MdcSegInfoAxialO *)aSeg->info())->phi0() < info->phi0()) {
00097 continue; }
00098 segList[isuper].insert(inSeg, iseg);
00099 isInserted = 1;
00100 break;
00101 }
00102 if (isInserted == 0) segList[isuper].append(inSeg);
00103 inSeg = (MdcSeg *) inSeg->next();
00104 }
00105
00106 }
00107
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119 }
00120
00121
00122 int MdcSegGrouperAx::incompWithSeg(const MdcSeg *refSeg,
00123 const MdcSeg *testSeg) {
00124
00125
00126
00127
00128 if (testSeg == 0) return 0;
00129 if(3 == _debug) {
00130
00131 testSeg->plotSegAll();
00132 }
00133
00134 MdcSegInfoAxialO *refInfo = (MdcSegInfoAxialO *) refSeg->info();
00135 MdcSegInfoAxialO *testInfo = (MdcSegInfoAxialO *) testSeg->info();
00136
00137
00138 double sigPhi0 = (refInfo->sigPhi0() > testInfo->sigPhi0() ?
00139 refInfo->sigPhi0() : testInfo->sigPhi0());
00140 double refPhi0 = refInfo->phi0();
00141 double testPhi0 = testInfo->phi0();
00142 double corrPhi0 = mdcWrapAng(refPhi0, testPhi0);
00143
00144 double sigCurv = (refInfo->sigCurv() > testInfo->sigCurv() ?
00145 refInfo->sigCurv() : testInfo->sigCurv());
00146 double refCurv = refInfo->curv();
00147 double testCurv = testInfo->curv();
00148
00149
00150 double phi0Cut = MdcTrkReconCut_combAxPhi0Cut;
00151 double curvCut = MdcTrkReconCut_combAxCurvCut;
00152
00153
00154
00155 if (g_tupleCombAx) {
00156 g_combAxdPhi0 = refPhi0 - corrPhi0;
00157 g_combAxdCurv = refCurv - testCurv;
00158 g_combAxQualitySeed = refSeg->quality()+refSeg->nHit()*10+refSeg->superlayer()->slayNum()*1000;
00159 g_combAxQualityTest = testSeg->quality()+testSeg->nHit()*10+testSeg->superlayer()->slayNum()*1000;
00160 g_combAxMcAmbigSeed = refSeg->testCombSegAmbig();
00161 g_combAxMcAmbigTest = testSeg->testCombSegAmbig();
00162 g_combAxMc = refSeg->testCombSeg(testSeg);
00163 g_combAxMcPt = refSeg->testCombSegPt();
00164 g_combAxMcTheta = refSeg->testCombSegTheta();
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178 g_tupleCombAx->write();
00179 }
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 if(refSeg->testCombSegPt()>0.4 && fabs(corrPhi0 - refPhi0)>0.4 &&(refSeg->testCombSeg(testSeg)>0.5)){
00193 if(3==_debug){
00194 std::cout<< endl<<" test " << std::endl;
00195 std::cout<<"seed seg: "; refSeg->plotSegAll(); std::cout<< std::endl;
00196 std::cout<<"test seg: "; testSeg->plotSegAll(); std::cout<< std::endl;
00197 std::cout<< " dPhi0 abnormal "<<corrPhi0 - refPhi0<<std::endl;
00198 }
00199 }else{
00200 if(3== _debug){
00201 std::cout<< endl<<" test " << std::endl;
00202 std::cout<<"seed seg: "; refSeg->plotSegAll(); std::cout<< std::endl;
00203 std::cout<<"test seg: "; testSeg->plotSegAll(); std::cout<< std::endl;
00204 std::cout<< " dPhi0 ok " <<setprecision(3)<< corrPhi0 - refPhi0<<std::endl;
00205 }
00206 }
00207
00208
00209
00210 cout<<setiosflags(ios::fixed);
00211
00212 if (fabs(corrPhi0 - refPhi0) > phi0Cut) {
00213 if(3 == _debug) std::cout << " SKIP by dPhi0 "
00214 <<fabs(corrPhi0-refPhi0)<<">"<<phi0Cut<<std::endl;
00215
00216
00217
00218 return -1;
00219 }else{
00220 if(3 == _debug)std::cout<< " dphi " <<setprecision(3)<< fabs(corrPhi0 - refPhi0);
00221 }
00222
00223
00224
00225
00226 if (fabs(refCurv - testCurv) > curvCut){
00227 if(3 == _debug) std::cout << " SKIP by dCurv"
00228 <<fabs(refCurv-testCurv)<<">"<<curvCut<<std::endl;
00229 return -2;
00230 }else{
00231 if(3 == _debug)std::cout<< " dCurv " <<setprecision(3)<< fabs(refCurv - testCurv);
00232 }
00233 if(3 == _debug) std::cout << " KEEP "<<std::endl;
00234
00235 cout<<setprecision(6);
00236 cout<<setiosflags(ios::scientific);
00237
00238 return 0;
00239 }
00240
00241
00242 int MdcSegGrouperAx::incompWithGroup(MdcSeg **segGroup, const MdcSeg *testSeg,
00243 int iply) {
00244
00245
00246 return 0;
00247 }
00248
00249
00250 void MdcSegGrouperAx::resetComb(const MdcSeg *seed) {
00251
00252
00253
00254 if (isValid != 0) {
00255 int i;
00256 for (i = 0; i < nDeep; i++) {
00257 delete [] isValid[i];
00258 isValid[i] = 0;
00259 }
00260 }
00261
00262 _seed = seed;
00263
00264 int islay = 0;
00265 int iply = 0;
00266 nPlyFilled = 0;
00267 nNull = 0;
00268 const MdcSuperLayer *seedSlay = 0;
00269 if (seed != 0) seedSlay = seed->superlayer();
00270
00271
00272
00273 for (const MdcSuperLayer *thisSlay = _gm->firstSlay(); thisSlay != 0;
00274 thisSlay = thisSlay->next()) {
00275 bool noGoodYet = true;
00276 islay++;
00277
00278 if (thisSlay == seedSlay) continue;
00279 if (thisSlay->whichView() != 0) continue;
00280 firstGood[iply] = 0;
00281
00282
00283 firstBad[iply] = 0;
00284 if (segList[islay-1].length() != 0) {
00285 isValid[iply] = new bool[segList[islay-1].length()];
00286 }
00287
00288 if(3 == _debug && segList[islay-1].length()>0) {
00289 std::cout<<std::endl<< "--match axial seg by phi in slayer "
00290 <<thisSlay->slayNum()<<"--"<<std::endl;
00291
00292
00293 }
00294
00295 for (int i = 0; i < (int) segList[islay-1].length(); i++) {
00296 MdcSeg *aSeg = segList[islay-1][i];
00297 int invalid = incompWithSeg(seed, aSeg);
00298 isValid[iply][i] = true;
00299 if (invalid < 0) {
00300 firstBad[iply] = i;
00301 isValid[iply][i] = false;
00302 if (noGoodYet) firstGood[iply] = i+1;
00303 } else if (invalid > 0) {
00304
00305 firstBad[iply] = i;
00306 for (int j = i; j < (int) segList[islay-1].length(); j++)
00307 isValid[iply][j] = false;
00308 break;
00309 } else {
00310 firstBad[iply] = i+1;
00311 noGoodYet = false;
00312 }
00313 }
00314
00315
00316 if (firstGood[iply] > firstBad[iply]) firstGood[iply] = firstBad[iply];
00317 if (firstGood[iply] == firstBad[iply]) {
00318
00319 delete [] isValid[iply];
00320 isValid[iply] = 0;
00321 continue;
00322 }
00323
00324 combList[iply] = &segList[islay-1];
00325 leaveGap[iply] = false;
00326 iply++;
00327 }
00328 nPlyFilled = iply;
00329 resetSegCounters();
00330 maxNull = nPlyFilled - 2;
00331 maxNull++;
00332 }
00333
00334 MdcTrack* MdcSegGrouperAx::storePar(MdcTrack* trk, double parms[2],
00335 double chisq, TrkContext& context,
00336 double t0) {
00337
00338
00339 assert(trk == 0);
00340 BesAngle foundPhi0(parms[0]);
00341
00342
00343
00344 TrkExchangePar par(0.0, foundPhi0.rad(), 2.*parms[1], 0.0, 0.0);
00345 return new MdcTrack(_gm->nSuper(), par, chisq, context, t0);
00346 }
00347
00348
00349
00350
00351
00352