105 {
106
107
108
109
110
111
112
113
114
115
117 std::cout << "=======Print before arbitrateHits=======" << std::endl;
118 }
119
120 int nDeleted = 0;
121 std::vector<MdcTrack*> trksToKill;
122 trksToKill.reserve(4);
123
125
126
127 int* usedInTrackNum =
new int [
nTrack()];
128
130
131 int *refitTrack =
new int [
nTrack()];
132 for (
int i = 0; i <
nTrack(); i++) {
133 refitTrack[i] = 0;
134 }
135
136
137 int itrack;
138 for (itrack = 0; itrack <
nTrack(); itrack++) {
140 if (atrack == 0) continue;
142 trkXRef[itrack] = atrack;
143 }
144
145 for (itrack = 0; itrack <
nTrack(); itrack++) {
146
147 if (8 ==
tkParam.
lPrint) std::cout<<
"arbitrate track No."<<itrack<< std::endl;
149 if (atrack == 0) continue;
151 int lRefit = 0;
152 int trackOld = -1;
154 assert (tkFit != 0);
156 assert (hitList != 0);
157restart:
158 for (
int ii = 0; ii <
nTrack(); ii++) usedInTrackNum[ii] = 0;
159
160
161 int nPrev = 0;
162 int nHitDeleted = 0;
163 int maxGapLength = 0;
164 int nGapGE2= 0;
165 int nGapGE3= 0;
166 int nHitInLayer[43];
167 int nDeleteInLayer[43];
168 for(int i=0;i<43;i++){
169 nHitInLayer[i]=0;
170 nDeleteInLayer[i]=0;
171 }
172 if(8 ==
tkParam.
lPrint) std::cout<<
"--arbitrate--"<<std::endl;
174 int nUsed = ihit->hit()->nUsedHits();
176 std::cout<<"nUsed="<<nUsed<<":";
177 ihit->hit()->printAll(std::cout);
178 }
180 double deltaChi = -999;
181 ihit->getFitStuff(deltaChi);
182 std::cout<< "deltaChi="<<deltaChi<<std::endl;
183 }
184 int layer = ihit->layerNumber();
185 nHitInLayer[layer]++;
186
187 if (!ihit->isActive()) {
188
189
190
192 nDeleteInLayer[layer]++;
194 std::cout<< "=remove above inactive "<<std::endl;
195 }
198 if(ihit == hitList->
end())
break;
199 --ihit;
200 }
201 continue;
202 }
203 if (nUsed > 1) {
204 bool wasUsed = false;
205 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator>
q =
206 ihit->hit()->getUsedHits();
208 if ( !i->isActive() ) continue;
210
211 int wasDel=0;
212 for(int idel = 0;idel<trksToKill.size();idel++){
213 if( recoTrk == &(trksToKill[idel]->track()) ) wasDel = 1;
214 }
215 if(wasDel==1) continue;
216
217 int id = recoTrk->
id();
218 if (
id == aRecoTrk.
id())
continue;
219 long index = 0;
220 bool findKey = idMap.
get(
id, index);
221
222 assert(index >= 0);
223 usedInTrackNum[index]++;
225 std::cout<<" track "<<itrack<<"&" <<index
226 << " shared hits "<<usedInTrackNum[index]<<":";
227 ihit->printAll(std::cout);
228 }
229 wasUsed = true;
230 }
231 if (wasUsed) nPrev++;
232 }
233 }
234
235 int testGap = 0;
236
237 for (int i=0;i<43;i++){
238
240 std::cout<<i<<" nHitInLayer "<<nHitInLayer[i]
241 <<" nDeleteInLayer "<<nDeleteInLayer[i]<<std::endl;
242 }
243
244 if(nHitInLayer[i]>0 && (nHitInLayer[i]-nDeleteInLayer[i])==0) {
245
246 nHitDeleted++;
248 cout << "rec hits have been deleted in this layer"<<std::endl;
249 }
250 testGap++;
251
252 }else if(nHitInLayer[i]==0){
253
254 testGap++;
255
256 }else{
257
258
259 if(testGap>=2){
260 nGapGE2++;
261 if(testGap>=3){ nGapGE3++; }
262 if(testGap>maxGapLength) maxGapLength=testGap;
263
264 }
265 testGap=0;
266 }
267 }
268
269 bool toBeDeleted = false;
270
271 if(
tkParam.
lPrint>1) std::cout<<
"arbitrateHits tkNo:"<<itrack<<
" nGapGE2= "<<nGapGE2 <<
" nGapGE3= "<<nGapGE3 <<
" maxGapLength= "<<maxGapLength<<std::endl;
272
273
277 <<" Killing tkNo " << itrack << endl;
278 }
279 toBeDeleted = true;
280 }
281
282
283
286 cout <<
"arbitrateHits: nGapGE2 "<<nGapGE2<<
" >= "<<
tkParam.
nGapGE2 <<
" Killing tkNo " << itrack << endl;
287 }
288 toBeDeleted = true;
289 }
292 cout <<
"arbitrateHits: nGapGE3 "<<nGapGE3<<
" >= "<<
tkParam.
nGapGE3 <<
" Killing tkNo " << itrack << endl;
293 }
294 toBeDeleted = true;
295 }
298 cout <<
"arbitrateHits: maxGapLength "<<maxGapLength<<
" >= "<<
tkParam.
maxGapLength<<
" Killing tkNo " << itrack << endl;
299 }
300 toBeDeleted = true;
301 }
302
303 if(toBeDeleted){
304 nDeleted++;
305 delete &(atrack->
track());
307 trksToKill.push_back(atrack);
308 continue;
309 }
310
311
312
313 int nMost = 0;
314 int trackMost = 0;
315 for (
int ii = 0; ii <
nTrack(); ii++) {
317 std::cout<<"tk:"<<itrack<<"&"<<ii
318 <<" shared "<<usedInTrackNum[ii]<<" hits "<< std::endl;
319 }
320 if (usedInTrackNum[ii] > nMost) {
321 nMost = usedInTrackNum[ii];
322 trackMost = ii;
323 }
324 }
325
326
327 if (trackMost == trackOld) {
328 std::cout << "ErrMsg(error) MdcTrackListBase:"
329 << "Something ghastly happened in MdcTrackListBase::arbitrateHits"
330 << std::endl;
331 return 0;
332 }
333 trackOld = trackMost;
334
335
336
337
338 double groupDiff = 0.0;
339
340 int nFound = 0;
343 int lGroupHits = 0;
344
347 std::cout<<"track "<<trackMost<<" shared "<<nMost<<" hits > Cut nOverlap "
349 }
350 lGroupHits = 1;
353 }
354
355
356
357
358
359 if(8 ==
tkParam.
lPrint) std::cout<<
"Go back through hits, looking up overlap hits"<< std::endl;
360 if (nMost > 0) {
363 int nUsed = ihit->hit()->nUsedHits();
364
366 std::cout<< "--hit go back, nUsed="<<nUsed<<":";
367 ihit->hit()->printAll(std::cout);
368 }
369
370
371 if (nUsed < 2) { continue; }
372
373
374 if (!ihit->isActive()) {
375 if (8 ==
tkParam.
lPrint){ std::cout<<
"act=0 continue"<<std::endl; }
376 continue;
377 }
378
379
380 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator>
q = ihit->hit()->getUsedHits();
381 while (
q.first!=
q.second) {
382 int dropThisHit = 0;
385
386 if (!otherHot->
isActive())
continue;
387
388
389 if ( &aRecoTrk == otherTrack) continue;
390 int otherId = otherTrack->
id();
391 long otherIndex = -1;
392 idMap.
get(otherId, otherIndex); assert(otherIndex >= 0);
393
394
395 if (lGroupHits && otherIndex != trackMost) continue;
396
397 if (lGroupHits) {
399 std::cout<<"group hits "<< std::endl;
400 }
401
402
403
404
405 int aDof = tkFit->
nActive() - 5;
408 if (aDof <= 0) {groupDiff = 999;}
409 else if (otherDof <= 0) {groupDiff = -999;}
410 else {
411 groupDiff += ihit->resid(0) * ihit->resid(0) * ihit->weight() /
412 aDof -
414 otherDof;
415 }
416 theseHits[nFound] =
const_cast<TrkHitOnTrk*
>(ihit.get());
417 thoseHits[nFound] = otherHot;
418 nFound++;
419 dropThisHit = 1;
420 } else {
421
423 std::cout<<"handle hits individually"<< std::endl;
424 }
425 nFound++;
426 if (fabs(ihit->resid(0)) > fabs(otherHot->
resid(0)) ) {
427
428 lRefit = 1;
429
430
431 const_cast<TrkHitOnTrk*
>(ihit.get())->setActivity(0);
432 dropThisHit = 1;
434 std::cout<<"dorp hit ";
435 const_cast<TrkHitOnTrk*
>(ihit.get())->print(std::cout);
436 }
437 break;
438 } else {
439
440 refitTrack[otherIndex] = 1;
441
444 std::cout<<"inactive hit on other track";
445 const_cast<TrkHitOnTrk*
>(ihit.get())->print(std::cout);
446 }
447 break;
448 }
449 }
450
451 if (dropThisHit == 1) break;
452
453 }
454
455
456 if (lGroupHits && nFound == nMost || nFound == nPrev) {
458 std::cout<<"we've found all of the shared hits on this track,Quit"<<std::endl;
459 }
460 break;
461 }
462
463 }
464
465
466 if (lGroupHits) {
468 cout << "nGroup: " << nMost << " groupDiff: " << groupDiff << endl;
469 cout <<
"Track: " << aRecoTrk.
id() <<
" nHit: "
470 << hitList->
nHit() <<
" nActive: "
471 << tkFit->
nActive() <<
" chisq/dof: " <<
474 cout <<
"Track: "<< othTrack.
id() <<
" nHit: " <<
475 othTrack.
hits()->
nHit() <<
" nActive: " <<
479 }
480
481 if (groupDiff > 0.0) {
482
483 lRefit = 1;
484 for (int ii = 0; ii < nMost; ii++) {
488
489 }
490 if (8 ==
tkParam.
lPrint) std::cout<<
"inactive hits on this track, No."<<aRecoTrk.
id()<< std::endl;
491 } else {
492
493 refitTrack[trackMost] = 1;
494 for (int ii = 0; ii < nMost; ii++) {
499
500
501 }
502 if (8 ==
tkParam.
lPrint) std::cout<<
"inactive hits on other track "<< std::endl;
503 }
504 delete [] theseHits;
505 delete [] thoseHits;
506
507 }
508
509 }
510
511
512
514 long index = -1;
515 idMap.
get(aRecoTrk.
id(), index); assert (index >= 0);
516
517 if (lRefit || refitTrack[index] == 1) {
519 std::cout<<
"after group ,refit track"<<aRecoTrk.
id()<< std::endl;
520 }
521 fitResult = hitList->
fit();
525 fitResult.
print(std::cerr);
526 }
527
528
529 double chisqperDOF;
530 bool badFit = true;
532 badFit = false;
533 int nDOF = tkFit->
nActive() - 5;
534 if (nDOF > 5){
535 chisqperDOF = tkFit->
chisq() / nDOF;
536 }else{
537 chisqperDOF = tkFit->
chisq();
538 }
539
542 double tem2 = (float) hitList->
nHit() - tkFit->
nActive();
546 badFit = true;
547 }
548 }
554 << std::endl;
555
556
557 }
559 cout <<
"Refitting track " << aRecoTrk.
id() <<
" success = "
560 << fitResult.
success() <<
"\n";
561 }
562
563 if (fitResult.
failure() || badFit ) {
564 nDeleted++;
565
566
567
569 cout <<
"fitResult.failure? "<<fitResult.
failure()
570 <<" badFit? "<<badFit <<" Killing tkNo " << itrack << endl;
571 }
572
573
574 trksToKill.push_back(atrack);
575 continue;
576 }
577 }
578
579 if (lGroupHits) goto restart;
580
581 }
582 if (8 ==
tkParam.
lPrint) std::cout<<
"end of loop over tracks"<< std::endl;
583
584
585 for (int itk = 0; itk < (int)trksToKill.size(); itk++) {
586 delete &(trksToKill[itk]->track());
587 trksToKill[itk]->setTrack(0);
589 if (8 ==
tkParam.
lPrint) std::cout<<
"remode dead track No."<<itk<< std::endl;
590 }
591 if (8 ==
tkParam.
lPrint) std::cout<<
"---end of arbitrateHits"<< std::endl;
592
593 delete [] usedInTrackNum;
594 delete [] refitTrack;
595 delete [] trkXRef;
596 return nDeleted;
597}
****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