296 {
297
298
299
300
301
302
303
304
305
308 cout << "Before pickHits";
309 if (is2d) cout<<" 2d:";
310 else cout<<" 3d:";
311 cout<< endl;
312 }
313
314 int nFound = 0;
315 int nCand = 0;
316 double rMin, rMax;
317 double rEnter, rExit;
319 int wireLow, wireHigh;
320 double phiLow, phiHigh;
321 int nCell;
324
325 double cellWidth;
326 double goodDriftCut;
327 double aresid = 0.;
328 int firstInputHit = -1;
329 int lCurl = 0;
330
331
334
337 assert (tkFit != 0);
339 assert (tkStatus != 0);
341 assert (hitList != 0);
343 double d0 = par.
d0();
344 double curv = 0.5 * par.
omega();
345 double sinPhi0 =
sin(par.
phi0());
346 double cosPhi0 =
cos(par.
phi0());
347
348
350 double absd0 = fabs(d0);
352
353 bool willCurl = false;
354 double rCurl = fabs(d0 + 1./curv);
356
357 if (rCurl < rMax) {
358 willCurl = true;
360 }
361
362 bool reachedLastInput = false;
363 bool hasCurled = false;
364
365
366 bool isHotOnLayer[43];
368 for(int ii=0; ii<43; ii++) isHotOnLayer[ii]=false;
370 isHotOnLayer[ihit->layerNumber()]=true;
371 }
372 }
373
374
375
377 layer = ((!lCurl) ? ( (hasCurled) ? gm->
prevLayer(layer) :
379
380
381 if (lCurl) {
382 lCurl = 0;
383 hasCurled = true;
384 }
385 if (tkStatus->
is2d() && layer->view() != 0)
continue;
386
388
389
390
391 bool lContinue = true;
393 if (layer == lastInputLayer) reachedLastInput = true;
394 }
395
396
397 if (hasCurled) {
398 rEnter = layer->rOut();
399 if (rEnter < rMin) {
400
401
402 break;
403 }
404 rExit = layer->rIn();
405 if (rExit < rMin) rExit = rMin;
406 if (rEnter > rMax) rEnter = rMax;
407
408
409
410 } else {
411 rEnter = layer->rIn();
412 rExit = layer->rOut();
413
414
415 if (rExit < rMin) continue;
416
417 if (willCurl) {
418 if (rEnter > rMax) {
419
420
421 hasCurled = 1;
422 continue;
423 }
424 if (rExit > rMax) {
425 lCurl = 1;
426 rExit = rMax;
427 }
428 } else {
429
430 if (rEnter > rMax) {
431
432 break;
433 }
434 if (rExit > rMax) rExit = rMax;
435 }
436 }
437
438 nCell = layer->nWires();
440
442
443 double deltaPhiCellWidth = 0.5 * (cellWidth * M_SQRT2)/layer->rMid();
444
445
447 int ierror = trk->
projectToR(rEnter, phiTemp, hasCurled);
448 phiEnter = phiTemp.
posRad();
449 if (ierror != 0) {
451 << "Error in MdcTrackList::pickHits projection, ierror " <<ierror<< std::endl;
452 continue;
453 }
454 ierror = trk->
projectToR(rExit, phiTemp, hasCurled);
455 phiExit = phiTemp.
posRad();
456 if (ierror != 0) {
458 << "Error in MdcTrackList::pickHits projection, ierror "<<ierror << std::endl;
459 continue;
460 }
461
462
464 std::cout<< endl<<"--pickHit of layer "<<layer->layNum()<<"--"<<std::endl;
465 std::cout<<
" track phiEnter "<< phiEnter.
deg()<<
" phiExit "<<phiExit.
deg()<<
" degree"<< std::endl;
468 std::cout<<
" goodDriftCut "<< goodDriftCut <<
"=sqrt(2)*0.5*"<<cellWidth<<
"+"<<
tkParam.
pickHitMargin<< std::endl;
469 }
470
472
473
474 double t_phiHit = -999.;
475 if (curv*Bz <= 0.0) {
476
477 phiLow = phiEnter;
478 phiHigh = phiExit;
479
482 } else {
483 phiLow = phiExit;
484 phiHigh = phiEnter;
487 }
488
489 if((phiHigh>0 && phiLow<0)){
491 }else if((phiHigh<0 && phiLow>0)){
493 }
494
495 double lowFloat = nCell /
Constants::twoPi * (phiLow - layer->phiOffset()) + 0.5;
496 double highFloat = nCell /
Constants::twoPi * (phiHigh - layer->phiOffset()) + 0.5;
497 wireLow = (lowFloat >= 0.0) ? int(lowFloat) : int(floor(lowFloat));
498 wireHigh = (highFloat >= 0.0) ? int(highFloat) : int(floor(highFloat));
499
501 std::cout <<
" wireLow "<<wireLow <<
" wireHigh "<<wireHigh <<
" phiLow "<<phiLow*180./
Constants::pi <<
" phiHigh "<<phiHigh*180./
Constants::pi << std::endl;
502 }
503
504 int tempDiff = 0;
507 if(wireLow>tempN/2. && wireHigh<tempN/2.){
509 tempDiff = wireHigh+tempN - wireLow +1;
510 }else{
512 tempDiff = wireHigh - wireLow +1;
513 }
514 }
515
516 if(wireLow>wireHigh) wireLow -= nCell;
517 long t_iHit = 0;
518 for (int thisWire = wireLow; thisWire <= wireHigh; thisWire++) {
520 thisHit = map->
hitWire(layer->layNum(), corrWire);
521
523 if(thisHit != 0) {cout<<endl<<
"test hit "; thisHit->
print(std::cout);}
524 else std::cout << endl<<"test hit ("<<layer->layNum()<<","<<corrWire<<")"<< std::endl;
525 }
526
527 double tof = 0.;
528 double z = 0.;
529 double driftDist = 0.;
530
531 double delx, dely;
532 double resid = 0., predDrift = 0.;
533 int ambig = 0;
535 double mcTkId = -9999;
536 if (thisHit == 0 ) {
537 delx = -d0 * sinPhi0 - layer->xWire(corrWire);
538 dely = d0 * cosPhi0 - layer->yWire(corrWire);
539 predDrift = cosPhi0 * dely - sinPhi0 * delx +
540 curv * (delx * delx + dely * dely);
541 if(6==
tkParam.
lPrint) cout<<
"No hit. predDrift="<<predDrift<<endl;
542 continue;
543 } else {
545
547 alink = 0;
548 }else{
551 }
552
553 if (alink == 0 || pickAm) {
554 if ((!tkStatus->
is2d()) && layer->view() != 0){
555 double rc = 9999.;
557 double rw = layer->rMid();
558 double alpha = acos(1 - rw*rw/(2*rc*rc));
559 double fltLen = rw;
560 if(fabs(1 - rw*rw/(2*rc*rc))<1) fltLen = rc *
alpha;
561 z = par.
z0() + fltLen* par.
tanDip();
563 double x = layer->getWire(corrWire)->xWireDC(z);
564 double y = layer->getWire(corrWire)->yWireDC(z);
565 delx = -d0 * sinPhi0 -
x;
566 dely = d0 * cosPhi0 - y;
568 }else{
569 delx = -d0 * sinPhi0 - thisHit->
x();
570 dely = d0 * cosPhi0 - thisHit->
y();
572 }
573 predDrift = cosPhi0 * dely - sinPhi0 * delx +
574 curv * (delx * delx + dely * dely);
575
576
577 ambig = (predDrift >= 0) ? 1 : -1;
578 if (hasCurled) ambig = -ambig;
579 double entranceAngle=0.;
580 driftDist = thisHit->
driftDist(tof+bunchTime,ambig,entranceAngle,0.,z);
581 if (alink == 0) {
582
583 resid = ambig * fabs(driftDist) - predDrift;
584 aresid = fabs(resid);
585
586 }
587 } else {
588 ambig = alink->
ambig();
591 }
592 }
593
594
595 double zGuess = par.
z0() + layer->rMid() * par.
tanDip();
596 double phiDCz = layer->getWire(corrWire)->phiDC(zGuess);
597 BesAngle phiDCzMax(phiDCz + deltaPhiCellWidth);
598 BesAngle phiDCzMin(phiDCz - deltaPhiCellWidth);
599
601 double sigma = 999.;
602 if (thisHit != 0 &&alink==0) {
603 double entranceAngle = 0.;
604 sigma = thisHit->
sigma(driftDist, ambig, entranceAngle, atan(par.
tanDip()), z);
605 }
609 if(curv*Bz<=0.){
610
611 if((phiDCzMin-phiExit>0) || (phiDCzMax-phiEnter<0))
m_pickPhizOk[t_iHit] = 0;
613 }else{
614 if((phiDCzMin-phiEnter>0) || (phiDCzMax-phiExit<0))
m_pickPhizOk[t_iHit] = 0;
616 }
622 double t_phiLowCut=-999.;
623 double t_phiHighCut= -999.;
624 if(t_phiHit > -998.){
625 t_phiLowCut = (phiEnter-t_phiHit)*rEnter;
626 t_phiHighCut = (phiExit-t_phiHit)*rExit;
627 }
634 t_iHit++;
635 }
636
637 if(curv*Bz<=0.){
638
639 if((phiDCzMin-phiExit>0) || (phiDCzMax-phiEnter<0)) {
640 if(6==
tkParam.
lPrint){ std::cout<<
" CUT by phiDCz not in phi En Ex range, curv>=0"<<std::endl; }
641 continue;
642 }
643 }else{
644
645 if((phiDCzMin-phiEnter>0) || (phiDCzMax-phiExit<0)) {
646 if(6==
tkParam.
lPrint){ std::cout<<
" CUT by phiDCz not in phi En Ex range, curv<0"<<std::endl; }
647 continue;
648 }
649 }
650
651
652
653 if (ambig != 0 && fabs(predDrift) > goodDriftCut){
654 if(6==
tkParam.
lPrint){cout<<
" predDrift "<<predDrift<<
">goodDriftCut "<<goodDriftCut<<endl;}
655 continue;
656 }
657
658
659 if (ambig == 0 && fabs(predDrift) > goodDriftCut){
661 cout<<" ambig==0 && |predDirft| "<<fabs(predDrift) <<"> goodDriftCut "<< goodDriftCut<<endl;
662 cout<<" No good hit, but track near cell-edge " << endl;
663 }
664 continue;
665 }
666
667
668
669 if (ambig != 0) {
670
671
672 double entranceAngle = 0.;
673 double sigma = thisHit->
sigma(driftDist, ambig, entranceAngle, atan(par.
tanDip()), z);
675
676
678
679
680
681
682 if (alink == 0 && (aresid <= residCut) ) {
684 cout << " (2) New hit found " << endl;
685 }
686
687 int isActive = 1;
688
689
690
691
692
694 tempHot.setActivity(isActive);
695
697 tempHot.setUsability(false);
698 if(6==
tkParam.
lPrint) std::cout<<
" thisHit used, setUsability false " << std::endl;
699 }
700
701 double flt = layer->rMid();
703 flt += 0.000001 * (thisHit->
x() + thisHit->
y());
704 tempHot.setFltLen(flt);
707 std::cout<< " Append Hot " << std::endl;
708 }
710 }else{
712 if(alink!=0){
713 std::cout<< "Exist hot found"<<std::endl;
714 }else if(aresid > residCut){
715 thisHit->
print(std::cout);
716 std::cout<<
" Drop hit. aresid "<<aresid<<
">residCut " <<residCut<<
" nSig "<<aresid/sigma<<
" nSigCut "<<(
tkParam.
maxActiveSigma*factor)<<
" factor "<<factor<<
" maxActiveSigma "<<
tkParam.
maxActiveSigma<<
" tanDip "<<par.
tanDip()<<std::endl;
717 }
718 }
719 }
720 if (!localHistory.member(
const_cast<MdcHitOnTrack*
>(alink))) {
723 nFound++;
724 if(6==
tkParam.
lPrint) std::cout<<
" nFound="<<nFound<<
" nCand "<<nCand<<std::endl;
725 if (layer == firstInputLayer && firstInputHit < 0) {
726 firstInputHit = nCand;
727 }
728 } else {
729 if(6==
tkParam.
lPrint) std::cout <<
"ErrMsg(warning) "<<
"would have inserted identical HOT "
730 "twice in history buffer" << std::endl;
731 }
732 }
733
734
735 else if (ambig == 0 && reachedLastInput) {
737 lContinue = false;
738 int nSuccess = 0;
739 int last = 8;
740 if (nCand < 8) last = nCand;
741 for (int i = 0; i < last; i++) {
742 int j = nCand - 1 - i;
743 if (localHistory[j] != 0) {
744
745 nSuccess++;
746 }
747 if (i == 2 && nSuccess >= 2) lContinue = true;
748 if (i == 4 && nSuccess >= 3) lContinue = true;
749 if (i == 7 && nSuccess >= 5) lContinue = true;
750 if(6==
tkParam.
lPrint) cout <<i<<
" (3) No hit found; if beyond known good region " << endl;
751 if (lContinue) {
752 if(6==
tkParam.
lPrint) std::cout<<
" pickHits break by lContinue. i="<<i<<
" nSuccess="<<nSuccess<< std::endl;
753 break;
754 }
755 }
756
757 if(6==
tkParam.
lPrint) cout <<
" (3) No hit found; if beyond known good region " << endl;
758
759 if (!lContinue) {
760
761 break;
762 }
763 localHistory.append(0);
764 }
765
766 nCand++;
767
768 if (ambig != 0 && reachedLastInput) {
772 }
773 }
774 else {
777 }
778 }
779 }
780
781 }
788 }
789 if (!lContinue && reachedLastInput) {
790
791 break;
792 }
793 }
794
795
796 bool lContinue = true;
797 for (int ihit = firstInputHit; ihit >= 0; ihit--) {
798 if (localHistory[ihit] != 0) {
799 if (lContinue) {
800
801 const MdcHitOnTrack *mdcHit = localHistory[ihit]->mdcHitOnTrack();
802 if (mdcHit != 0) {
805 }
806 }
807 continue;
808 } else {
809
811 std::cout << " gap found; delete hits. ";
812 }
813 if (!localHistory[ihit]->isUsable()) {
815 cout << "about to delete hit for unusable HOT:" << endl;
816 localHistory[ihit]->print(std::cout);
817 }
818 hitList->
removeHit(localHistory[ihit]->hit());
819 }
821 cout << " current contents of localHistory: "
822 <<localHistory.length()<<"hot" << endl;
823
824
825
826 }
827 nFound--;
828 }
829 }
830 else if (localHistory[ihit] == 0) {
831 if(6==
tkParam.
lPrint){ cout <<
" localHistory= 0 " << endl; }
832 int nSuccess = 0;
833 lContinue = false;
834 int last = 8;
835 if (nCand < 8) last = nCand;
836 for (int i = 0; i < last; i++) {
837 int j = ihit + 1 + i;
838 if (localHistory[j] != 0) nSuccess++;
839 if (i == 2 && nSuccess >= 2) lContinue = true;
840 if (i == 4 && nSuccess >= 3) lContinue = true;
841
842 if (lContinue) break;
843 }
844 }
845 }
848 }
850 if(6==
tkParam.
lPrint){ cout <<
" localHistory= 0 " << endl; }
851
853 cout << "After pickHits found " << nFound <<" hit(s)"<< endl;
855 std::cout<< std::endl;
856 }
857
858 return nFound;
859}
double sin(const BesAngle a)
double cos(const BesAngle a)
double abs(const EvtComplex &c)
NTuple::Array< long > m_pickIs2D
NTuple::Array< double > m_pickPt
NTuple::Item< long > m_pickLayer
NTuple::Array< long > m_pickMcTk
AIDA::IHistogram1D * g_pickHitWire
NTuple::Item< long > m_pickNCellWithDigi
NTuple::Array< double > m_pickPhiLowCut
double haveDigiDrift[43][288]
NTuple::Array< double > m_pickDriftCut
NTuple::Array< double > m_pickCurv
NTuple::Array< double > m_pickDrift
NTuple::Array< double > m_pickResid
NTuple::Array< int > m_pickPhizOk
NTuple::Array< double > m_pickSigma
NTuple::Array< long > m_pickWire
NTuple::Tuple * m_tuplePick
NTuple::Array< double > m_pickDriftTruth
NTuple::Array< double > m_pickZ
NTuple::Array< double > m_pickPredDrift
NTuple::Item< long > m_pickNCell
NTuple::Array< double > m_pickPhiHighCut
static const double twoPi
static const int maxCell[43]
static const double radToDegrees
const MdcLayer * prevLayer(int lay) const
const MdcLayer * lastLayer() const
const MdcLayer * firstLayer() const
const MdcLayer * nextLayer(int lay) const
MdcHit * hitWire(int lay, int wire) const
const MdcLayer * layer() const
const MdcDigi * digi() const
void print(std::ostream &o) const
double sigma(double, int, double, double, double) const
double driftDist(double, int, double, double, double) const
const MdcLayer * layer() const
void setFirstLayer(const MdcLayer *l)
int projectToR(double radius, BesAngle &phiIntersect, int lCurl=0) const
void setHasCurled(bool c=true)
const MdcLayer * firstLayer() const
const MdcLayer * lastLayer() const
void setLastLayer(const MdcLayer *l)
int getTrackIndex() const
const TrkHitOnTrk * getHitOnTrack(const TrkRecoTrk *trk) const
hot_iterator begin() const
TrkHitOnTrk * appendHot(const TrkHitOnTrk *theHot)
bool removeHit(const TrkFundHit *theHit)
void print(std::ostream &o) const
const BField & bField() const
int mdcWrapWire(int wireIn, int nCell)