104 {
105
106
107
108
109
110
111
112
113
114
116 std::cout << "=======Print before arbitrateHits=======" << std::endl;
117 }
118
119 int nDeleted = 0;
120 std::vector<MdcTrack*> trksToKill;
121 trksToKill.reserve(4);
122
124
125
126 int* usedInTrackNum =
new int [
nTrack()];
127
129
130 int *refitTrack =
new int [
nTrack()];
131 for (
int i = 0; i <
nTrack(); i++) {
132 refitTrack[i] = 0;
133 }
134
135
136 int itrack;
137 for (itrack = 0; itrack <
nTrack(); itrack++) {
139 if (atrack == 0) continue;
141 trkXRef[itrack] = atrack;
142 }
143
144 for (itrack = 0; itrack <
nTrack(); itrack++) {
145
146 if (8 ==
tkParam.
lPrint) std::cout<<
"arbitrate track No."<<itrack<< std::endl;
148 if (atrack == 0) continue;
150 int lRefit = 0;
151 int trackOld = -1;
153 assert (tkFit != 0);
155 assert (hitList != 0);
156restart:
157 for (
int ii = 0; ii <
nTrack(); ii++) usedInTrackNum[ii] = 0;
158
159
160 int nPrev = 0;
161 int nHitDeleted = 0;
162 int maxGapLength = 0;
163 int nGapGE2= 0;
164 int nGapGE3= 0;
165 int nHitInLayer[43];
166 int nDeleteInLayer[43];
167 for(int i=0;i<43;i++){
168 nHitInLayer[i]=0;
169 nDeleteInLayer[i]=0;
170 }
171 if(8 ==
tkParam.
lPrint) std::cout<<
"--arbitrate--"<<std::endl;
173 int nUsed = ihit->hit()->nUsedHits();
175 std::cout<<"nUsed="<<nUsed<<":";
176 ihit->hit()->printAll(std::cout);
177 }
179 double deltaChi = -999;
180 ihit->getFitStuff(deltaChi);
181 std::cout<< "deltaChi="<<deltaChi<<std::endl;
182 }
183 int layer = ihit->layerNumber();
184 nHitInLayer[layer]++;
185
186 if (!ihit->isActive()) {
187
188
189
191 nDeleteInLayer[layer]++;
193 std::cout<< "=remove above inactive "<<std::endl;
194 }
197 if(ihit == hitList->
end())
break;
198 --ihit;
199 }
200 continue;
201 }
202 if (nUsed > 1) {
203 bool wasUsed = false;
204 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator>
q =
205 ihit->hit()->getUsedHits();
207 if ( !i->isActive() ) continue;
209 int id = recoTrk->
id();
210 if (
id == aRecoTrk.
id())
continue;
211 long index = 0;
212 idMap.
get(
id, index);
213 assert(index >= 0);
214 usedInTrackNum[index]++;
216 std::cout<<" track "<<itrack<<"&" <<index
217 << " shared hits "<<usedInTrackNum[index]<<":";
218 ihit->printAll(std::cout);
219 }
220 wasUsed = true;
221 }
222 if (wasUsed) nPrev++;
223 }
224 }
225
226 int testGap = 0;
227
228 for (int i=0;i<43;i++){
229
231 std::cout<<i<<" nHitInLayer "<<nHitInLayer[i]
232 <<" nDeleteInLayer "<<nDeleteInLayer[i]<<std::endl;
233 }
234
235 if(nHitInLayer[i]>0 && (nHitInLayer[i]-nDeleteInLayer[i])==0) {
236
237 nHitDeleted++;
239 cout << "rec hits have been deleted in this layer"<<std::endl;
240 }
241 testGap++;
242
243 }else if(nHitInLayer[i]==0){
244
245 testGap++;
246
247 }else{
248
249
250 if(testGap>=2){
251 nGapGE2++;
252 if(testGap>=3){ nGapGE3++; }
253 if(testGap>maxGapLength) maxGapLength=testGap;
254
255 }
256 testGap=0;
257 }
258 }
259
260 bool toBeDeleted = false;
261
262 if(
tkParam.
lPrint>1) std::cout<<
"arbitrateHits tkNo:"<<itrack<<
" nGapGE2= "<<nGapGE2 <<
" nGapGE3= "<<nGapGE3 <<
" maxGapLength= "<<maxGapLength<<std::endl;
263
264
268 <<" Killing tkNo " << itrack << endl;
269 }
270 toBeDeleted = true;
271 }
272
273
274
277 cout <<
"arbitrateHits: nGapGE2 "<<nGapGE2<<
" >= "<<
tkParam.
nGapGE2 <<
" Killing tkNo " << itrack << endl;
278 }
279 toBeDeleted = true;
280 }
283 cout <<
"arbitrateHits: nGapGE3 "<<nGapGE3<<
" >= "<<
tkParam.
nGapGE3 <<
" Killing tkNo " << itrack << endl;
284 }
285 toBeDeleted = true;
286 }
289 cout <<
"arbitrateHits: maxGapLength "<<maxGapLength<<
" >= "<<
tkParam.
maxGapLength<<
" Killing tkNo " << itrack << endl;
290 }
291 toBeDeleted = true;
292 }
293
294 if(toBeDeleted){
295 nDeleted++;
296 delete &(atrack->
track());
298 trksToKill.push_back(atrack);
299 continue;
300 }
301
302
303
304 int nMost = 0;
305 int trackMost = 0;
306 for (
int ii = 0; ii <
nTrack(); ii++) {
308 std::cout<<"tk:"<<itrack<<"&"<<ii
309 <<" shared "<<usedInTrackNum[ii]<<" hits "<< std::endl;
310 }
311 if (usedInTrackNum[ii] > nMost) {
312 nMost = usedInTrackNum[ii];
313 trackMost = ii;
314 }
315 }
316
317
318 if (trackMost == trackOld) {
319 std::cout << "ErrMsg(error) MdcTrackListBase:"
320 << "Something ghastly happened in MdcTrackListBase::arbitrateHits"
321 << std::endl;
322 return 0;
323 }
324 trackOld = trackMost;
325
326
327
328
329 double groupDiff = 0.0;
330
331 int nFound = 0;
334 int lGroupHits = 0;
335
338 std::cout<<"track "<<trackMost<<" shared "<<nMost<<" hits > Cut nOverlap "
340 }
341 lGroupHits = 1;
344 }
345
346
347
348
349
350 if(8 ==
tkParam.
lPrint) std::cout<<
"Go back through hits, looking up overlap hits"<< std::endl;
351 if (nMost > 0) {
354 int nUsed = ihit->hit()->nUsedHits();
355
357 std::cout<< "--hit go back, nUsed="<<nUsed<<":";
358 ihit->hit()->printAll(std::cout);
359 }
360
361
362 if (nUsed < 2) { continue; }
363
364
365 if (!ihit->isActive()) {
366 if (8 ==
tkParam.
lPrint){ std::cout<<
"act=0 continue"<<std::endl; }
367 continue;
368 }
369
370
371 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator>
q = ihit->hit()->getUsedHits();
372 while (
q.first!=
q.second) {
373 int dropThisHit = 0;
376
377 if (!otherHot->
isActive())
continue;
378
379
380 if ( &aRecoTrk == otherTrack) continue;
381 int otherId = otherTrack->
id();
382 long otherIndex = -1;
383 idMap.
get(otherId, otherIndex); assert(otherIndex >= 0);
384
385
386 if (lGroupHits && otherIndex != trackMost) continue;
387
388 if (lGroupHits) {
390 std::cout<<"group hits "<< std::endl;
391 }
392
393
394
395
396 int aDof = tkFit->
nActive() - 5;
399 if (aDof <= 0) {groupDiff = 999;}
400 else if (otherDof <= 0) {groupDiff = -999;}
401 else {
402 groupDiff += ihit->resid(0) * ihit->resid(0) * ihit->weight() /
403 aDof -
405 otherDof;
406 }
407 theseHits[nFound] =
const_cast<TrkHitOnTrk*
>(ihit.get());
408 thoseHits[nFound] = otherHot;
409 nFound++;
410 dropThisHit = 1;
411 } else {
412
414 std::cout<<"handle hits individually"<< std::endl;
415 }
416 nFound++;
417 if (fabs(ihit->resid(0)) > fabs(otherHot->
resid(0)) ) {
418
419 lRefit = 1;
420
421
422 const_cast<TrkHitOnTrk*
>(ihit.get())->setActivity(0);
423 dropThisHit = 1;
425 std::cout<<"dorp hit ";
426 const_cast<TrkHitOnTrk*
>(ihit.get())->print(std::cout);
427 }
428 break;
429 } else {
430
431 refitTrack[otherIndex] = 1;
432
435 std::cout<<"inactive hit on other track";
436 const_cast<TrkHitOnTrk*
>(ihit.get())->print(std::cout);
437 }
438 break;
439 }
440 }
441
442 if (dropThisHit == 1) break;
443
444 }
445
446
447 if (lGroupHits && nFound == nMost || nFound == nPrev) {
449 std::cout<<"we've found all of the shared hits on this track,Quit"<<std::endl;
450 }
451 break;
452 }
453
454 }
455
456
457 if (lGroupHits) {
459 cout << "nGroup: " << nMost << " groupDiff: " << groupDiff << endl;
460 cout <<
"Track: " << aRecoTrk.
id() <<
" nHit: "
461 << hitList->
nHit() <<
" nActive: "
462 << tkFit->
nActive() <<
" chisq/dof: " <<
465 cout <<
"Track: "<< othTrack.
id() <<
" nHit: " <<
466 othTrack.
hits()->
nHit() <<
" nActive: " <<
470 }
471
472 if (groupDiff > 0.0) {
473
474 lRefit = 1;
475 for (int ii = 0; ii < nMost; ii++) {
479
480 }
481 if (8 ==
tkParam.
lPrint) std::cout<<
"inactive hits on this track, No."<<aRecoTrk.
id()<< std::endl;
482 } else {
483
484 refitTrack[trackMost] = 1;
485 for (int ii = 0; ii < nMost; ii++) {
489
490 }
491 if (8 ==
tkParam.
lPrint) std::cout<<
"inactive hits on other track "<< std::endl;
492 }
493 delete [] theseHits;
494 delete [] thoseHits;
495
496 }
497
498 }
499
500
501
503 long index = -1;
504 idMap.
get(aRecoTrk.
id(), index); assert (index >= 0);
505
506 if (lRefit || refitTrack[index] == 1) {
508 std::cout<<
"after group ,refit track"<<aRecoTrk.
id()<< std::endl;
509 }
510 fitResult = hitList->
fit();
514 fitResult.
print(std::cerr);
515 }
516
517
518 double chisqperDOF;
519 bool badFit = true;
521 badFit = false;
522 int nDOF = tkFit->
nActive() - 5;
523 if (nDOF > 5){
524 chisqperDOF = tkFit->
chisq() / nDOF;
525 }else{
526 chisqperDOF = tkFit->
chisq();
527 }
528
531 double tem2 = (float) hitList->
nHit() - tkFit->
nActive();
535 badFit = true;
536 }
537 }
543 << std::endl;
544
545
546 }
548 cout <<
"Refitting track " << aRecoTrk.
id() <<
" success = "
549 << fitResult.
success() <<
"\n";
550 }
551
552 if (fitResult.
failure() || badFit ) {
553 nDeleted++;
554
555
556
558 cout <<
"fitResult.failure? "<<fitResult.
failure()
559 <<" badFit? "<<badFit <<" Killing tkNo " << itrack << endl;
560 }
561 delete &(atrack->
track());
563 trksToKill.push_back(atrack);
564 continue;
565 }
566 }
567
568 if (lGroupHits) goto restart;
569
570 }
571 if (8 ==
tkParam.
lPrint) std::cout<<
"end of loop over tracks"<< std::endl;
572
573
574 for (int itk = 0; itk < (int)trksToKill.size(); itk++) {
576 if (8 ==
tkParam.
lPrint) std::cout<<
"remode dead track No."<<itk<< std::endl;
577 }
578 if (8 ==
tkParam.
lPrint) std::cout<<
"---end of arbitrateHits"<< std::endl;
579
580 delete [] usedInTrackNum;
581 delete [] refitTrack;
582 delete [] trkXRef;
583 return nDeleted;
584}
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
bool get(const K &theKey, V &theAnswer) const
void put(const K &, const V &)
void remove(MdcTrack *atrack)
void setTrack(TrkRecoTrk *trk)
virtual double chisq() const =0
void print(std::ostream &ostr) const
virtual void addHistory(const TrkErrCode &status, const char *modulename)
virtual int nActive() const =0
hot_iterator begin() const
bool removeHit(const TrkFundHit *theHit)
void setActivity(bool turnOn)
double resid(bool exclude=false) const
TrkRecoTrk * parentTrack() const
const TrkFundHit * hit() const
const TrkFit * fitResult() const
const TrkFitStatus * status() const