108 {
109
110
111
112
113
114
115
116
117
118
120 std::cout << "=======Print before arbitrateHits=======" << std::endl;
121 }
122
123 int nDeleted = 0;
124 std::vector<MdcTrack*> trksToKill;
125 trksToKill.reserve(4);
126
128
129
130 int* usedInTrackNum =
new int [
nTrack()];
131
133
134 int *refitTrack =
new int [
nTrack()];
135 for (
int i = 0; i <
nTrack(); i++) {
136 refitTrack[i] = 0;
137 }
138
139
140 int itrack;
141 for (itrack = 0; itrack <
nTrack(); itrack++) {
143 if (atrack == 0) continue;
145 trkXRef[itrack] = atrack;
146 }
147
148 for (itrack = 0; itrack <
nTrack(); itrack++) {
149
150 if (8 ==
tkParam.
lPrint) std::cout<<
"arbitrate track No."<<itrack<< std::endl;
152 if (atrack == 0) continue;
154 int lRefit = 0;
155 int trackOld = -1;
157 assert (tkFit != 0);
159 assert (hitList != 0);
160restart:
161 for (
int ii = 0; ii <
nTrack(); ii++) usedInTrackNum[ii] = 0;
162
163
164 int nPrev = 0;
165 int nHitDeleted = 0;
166 int maxGapLength = 0;
167 int nGapGE2= 0;
168 int nGapGE3= 0;
169 int nHitInLayer[43];
170 int nDeleteInLayer[43];
171 for(int i=0;i<43;i++){
172 nHitInLayer[i]=0;
173 nDeleteInLayer[i]=0;
174 }
175 if(8 ==
tkParam.
lPrint) std::cout<<
"--arbitrate--"<<std::endl;
177
178
180 int nUsed = ihit->hit()->nUsedHits();
182 std::cout<<"nUsed="<<nUsed<<":";
183 ihit->hit()->printAll(std::cout);
184 }
186 double deltaChi = -999;
187 ihit->getFitStuff(deltaChi);
188 std::cout<< "deltaChi="<<deltaChi<<std::endl;
189 }
190 int layer = ihit->layerNumber();
191 nHitInLayer[layer]++;
192
193 if (!ihit->isActive()) {
194
195
196
198 nDeleteInLayer[layer]++;
200 std::cout<< "=remove above inactive "<<std::endl;
201 }
204 if(ihit == hitList->
end())
break;
205 --ihit;
206 }
207 continue;
208 }
209 if (nUsed > 1) {
210 bool wasUsed = false;
211 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator>
q =
212 ihit->hit()->getUsedHits();
214 if ( !i->isActive() ) continue;
216
217 int wasDel=0;
218 for(int idel = 0;idel<trksToKill.size();idel++){
219 if( recoTrk == &(trksToKill[idel]->track()) ) wasDel = 1;
220 }
221 if(wasDel==1) continue;
222
223 int id = recoTrk->
id();
224 if (
id == aRecoTrk.
id())
continue;
225 long index = 0;
226 bool findKey = idMap.
get(
id, index);
227
228 assert(index >= 0);
229 usedInTrackNum[index]++;
231 std::cout<<" track "<<itrack<<"&" <<index
232 << " shared hits "<<usedInTrackNum[index]<<":";
233 ihit->printAll(std::cout);
234 }
235 wasUsed = true;
236 }
237 if (wasUsed) nPrev++;
238 }
239 }
240
241 int testGap = 0;
242 int first_layer = 0;
244
245 for (int i=first_layer;i<43;i++){
246
248 std::cout<<i<<" nHitInLayer "<<nHitInLayer[i]
249 <<" nDeleteInLayer "<<nDeleteInLayer[i]<<std::endl;
250 }
251
252 if(nHitInLayer[i]>0 && (nHitInLayer[i]-nDeleteInLayer[i])==0) {
253
254 nHitDeleted++;
256 cout << "rec hits have been deleted in this layer"<<std::endl;
257 }
258 testGap++;
259
260 }else if(nHitInLayer[i]==0){
261
262 testGap++;
263
264 }else{
265
266
267 if(testGap>=2){
268 nGapGE2++;
269 if(testGap>=3){ nGapGE3++; }
270 if(testGap>maxGapLength) maxGapLength=testGap;
271
272 }
273 testGap=0;
274 }
275 }
276
277 bool toBeDeleted = false;
278
279 if(
tkParam.
lPrint>1) std::cout<<
"arbitrateHits tkNo:"<<itrack<<
" nGapGE2= "<<nGapGE2 <<
" nGapGE3= "<<nGapGE3 <<
" maxGapLength= "<<maxGapLength<<std::endl;
280
281
285 <<" Killing tkNo " << itrack << endl;
286 }
287 toBeDeleted = true;
288 }
289
290
291
294 cout <<
"arbitrateHits: nGapGE2 "<<nGapGE2<<
" >= "<<
tkParam.
nGapGE2 <<
" Killing tkNo " << itrack << endl;
295 }
296 toBeDeleted = true;
297 }
300 cout <<
"arbitrateHits: nGapGE3 "<<nGapGE3<<
" >= "<<
tkParam.
nGapGE3 <<
" Killing tkNo " << itrack << endl;
301 }
302 toBeDeleted = true;
303 }
306 cout <<
"arbitrateHits: maxGapLength "<<maxGapLength<<
" >= "<<
tkParam.
maxGapLength<<
" Killing tkNo " << itrack << endl;
307 }
308 toBeDeleted = true;
309 }
310
311 if(toBeDeleted){
312 nDeleted++;
313 delete &(atrack->
track());
315 trksToKill.push_back(atrack);
316 continue;
317 }
318
319
320
321 int nMost = 0;
322 int trackMost = 0;
323 for (
int ii = 0; ii <
nTrack(); ii++) {
325 std::cout<<"tk:"<<itrack<<"&"<<ii
326 <<" shared "<<usedInTrackNum[ii]<<" hits "<< std::endl;
327 }
328 if (usedInTrackNum[ii] > nMost) {
329 nMost = usedInTrackNum[ii];
330 trackMost = ii;
331 }
332 }
333
334
335 if (trackMost == trackOld) {
336 std::cout << "ErrMsg(error) MdcTrackListBase:"
337 << "Something ghastly happened in MdcTrackListBase::arbitrateHits"
338 << std::endl;
339 return 0;
340 }
341 trackOld = trackMost;
342
343
344
345
346 double groupDiff = 0.0;
347
348 int nFound = 0;
351 int lGroupHits = 0;
352
355 std::cout<<"track "<<trackMost<<" shared "<<nMost<<" hits > Cut nOverlap "
357 }
358 lGroupHits = 1;
361 }
362
363
364
365
366
367 if(8 ==
tkParam.
lPrint) std::cout<<
"Go back through hits, looking up overlap hits"<< std::endl;
368 if (nMost > 0) {
372 int nUsed = ihit->hit()->nUsedHits();
373
375 std::cout<< "--hit go back, nUsed="<<nUsed<<":";
376 ihit->hit()->printAll(std::cout);
377 }
378
379
380 if (nUsed < 2) { continue; }
381
382
383 if (!ihit->isActive()) {
384 if (8 ==
tkParam.
lPrint){ std::cout<<
"act=0 continue"<<std::endl; }
385 continue;
386 }
387
388
389 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator>
q = ihit->hit()->getUsedHits();
390 while (
q.first!=
q.second) {
391 int dropThisHit = 0;
394
395 if (!otherHot->
isActive())
continue;
396
397
398 if ( &aRecoTrk == otherTrack) continue;
399 int otherId = otherTrack->
id();
400 long otherIndex = -1;
401 idMap.
get(otherId, otherIndex); assert(otherIndex >= 0);
402
403
404 if (lGroupHits && otherIndex != trackMost) continue;
405
406 if (lGroupHits) {
408 std::cout<<"group hits "<< std::endl;
409 }
410
411
412
413
414 int aDof = tkFit->
nActive() - 5;
417 if (aDof <= 0) {groupDiff = 999;}
418 else if (otherDof <= 0) {groupDiff = -999;}
419 else {
420 groupDiff += ihit->resid(0) * ihit->resid(0) * ihit->weight() /
421 aDof -
423 otherDof;
424 }
425 theseHits[nFound] =
const_cast<TrkHitOnTrk*
>(ihit.get());
426 thoseHits[nFound] = otherHot;
427 nFound++;
428 dropThisHit = 1;
429 } else {
430
432 std::cout<<"handle hits individually"<< std::endl;
433 }
434 nFound++;
435 if (fabs(ihit->resid(0)) > fabs(otherHot->
resid(0)) ) {
436
437 lRefit = 1;
438
439
441 dropThisHit = 1;
443 std::cout<<"dorp hit ";
445 }
446 break;
447 } else {
448
449 refitTrack[otherIndex] = 1;
450
453 std::cout<<"inactive hit on other track";
455 }
456 break;
457 }
458 }
459
460 if (dropThisHit == 1) break;
461
462 }
463
464
465 if (lGroupHits && nFound == nMost || nFound == nPrev) {
467 std::cout<<"we've found all of the shared hits on this track,Quit"<<std::endl;
468 }
469 break;
470 }
471
472 }
473
474
475 if (lGroupHits) {
477 cout << "nGroup: " << nMost << " groupDiff: " << groupDiff << endl;
478 cout <<
"Track: " << aRecoTrk.
id() <<
" nHit: "
479 << hitList->
nHit() <<
" nActive: "
480 << tkFit->
nActive() <<
" chisq/dof: " <<
483 cout <<
"Track: "<< othTrack.
id() <<
" nHit: " <<
484 othTrack.
hits()->
nHit() <<
" nActive: " <<
487 (othTrack.
fitResult()->nActive() - 5) << endl;
488 }
489
490 if (groupDiff > 0.0) {
491
492 lRefit = 1;
493 for (int ii = 0; ii < nMost; ii++) {
497
498 }
499 if (8 ==
tkParam.
lPrint) std::cout<<
"inactive hits on this track, No."<<aRecoTrk.
id()<< std::endl;
500 } else {
501
502 refitTrack[trackMost] = 1;
503 for (int ii = 0; ii < nMost; ii++) {
508
509
510 }
511 if (8 ==
tkParam.
lPrint) std::cout<<
"inactive hits on other track "<< std::endl;
512 }
513 delete [] theseHits;
514 delete [] thoseHits;
515
516 }
517
518 }
519
520
521
523 long index = -1;
524 idMap.
get(aRecoTrk.
id(), index); assert (index >= 0);
525
526 if (lRefit || refitTrack[index] == 1) {
528 std::cout<<
"after group ,refit track"<<aRecoTrk.
id()<< std::endl;
529 }
530 fitResult = hitList->
fit();
534 fitResult.
print(std::cerr);
535 }
536
537
538 double chisqperDOF;
539 bool badFit = true;
541 badFit = false;
542 int nDOF = tkFit->
nActive() - 5;
543 if (nDOF > 5){
544 chisqperDOF = tkFit->
chisq() / nDOF;
545 }else{
546 chisqperDOF = tkFit->
chisq();
547 }
548
551 double tem2 = (float) hitList->
nHit() - tkFit->
nActive();
555 badFit = true;
556 }
557 }
563 << std::endl;
564
565
566 }
568 cout <<
"Refitting track " << aRecoTrk.
id() <<
" success = "
569 << fitResult.
success() <<
"\n";
570 }
571
572 if (fitResult.
failure() || badFit ) {
573 nDeleted++;
574
575
576
578 cout <<
"fitResult.failure? "<<fitResult.
failure()
579 <<" badFit? "<<badFit <<" Killing tkNo " << itrack << endl;
580 }
581
582
583 trksToKill.push_back(atrack);
584 continue;
585 }
586 }
587
588 if (lGroupHits) goto restart;
589
590 }
591 if (8 ==
tkParam.
lPrint) std::cout<<
"end of loop over tracks"<< std::endl;
592
593
594 for (int itk = 0; itk < (int)trksToKill.size(); itk++) {
595 delete &(trksToKill[itk]->track());
596 trksToKill[itk]->setTrack(0);
598 if (8 ==
tkParam.
lPrint) std::cout<<
"remode dead track No."<<itk<< std::endl;
599 }
600 if (8 ==
tkParam.
lPrint) std::cout<<
"---end of arbitrateHits"<< std::endl;
601
602 delete [] usedInTrackNum;
603 delete [] refitTrack;
604 delete [] trkXRef;
605 return nDeleted;
606}
****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
virtual void print(std::ostream &) const
const TrkFit * fitResult() const
const TrkFitStatus * status() const