23 #include "Utilities/AssociationUtil.h" 24 #include "NovaDAQConventions/DAQConventions.h" 79 std::vector< std::vector<double> > &iP,
80 std::vector< bool> &propDir,
81 std::vector< std::vector<bool> > &trkHits,
82 std::vector< std::vector<std::array<double, 2> > > &trkTraj,
83 std::vector< std::vector<std::array<double, 2> > > &trkDirs );
87 std::vector< bool> &propDir,
88 std::vector< std::vector<bool> > &trkHits,
89 std::vector< std::vector<std::array<double, 2> > > &trkTraj,
90 std::vector< std::vector<std::array<double, 2> > > &trkDirs,
96 int zDir,
bool zProp, std::vector<bool> &trkHits,
97 std::vector<std::array<double, 2> > &trkTrajs,
98 std::vector<std::array<double, 2> > &trkDirs,
99 std::vector<double> &iP,
int &outnhits,
int firsthit = -1);
104 std::vector< std::vector<std::array<double, 2> > > &trkTrajs,
105 std::vector< std::vector<std::array<double, 2> > >&trkDirs,
106 std::vector< double > &
chi2,
107 std::vector< std::vector<double> > &pfiltState,
108 std::vector<bool> &zProp,
109 std::vector<int> &ntrkHits,
110 std::vector<bool> &ignoreTrk);
118 double FilterOnly(
const std::vector<bool> & trkHits,
119 std::vector<std::array<double, 2> > &trkTrajs,
120 std::vector<std::array<double, 2> > &trkDirs,
121 std::vector<double> &initialP,
int zDir,
bool zProp,
127 std::vector<bool> &trkHits,
128 std::vector<std::array<double, 2> > trkTrajs,
129 bool zProp, std::vector<bool> &newTrkHits);
133 void PlaneToCellSort(std::vector<std::array<double, 2> > &planeSortedHits);
134 void CellToPlaneSort(std::vector<std::array<double, 2> > &cellSortedHits);
136 int CountHits(
const std::vector<bool> & trkHits);
138 int GetFirstHit(
const std::vector<bool> & trkHits);
139 int GetLastHit(
const std::vector<bool> & trkHits);
150 produces< std::vector<rb::Track> >();
151 produces< art::Assns<rb::Track, rb::Cluster> >();
185 std::unique_ptr< std::vector<rb::Track> > trackcol(
new std::vector<rb::Track>);
196 for(
unsigned int i = 0;
i < slicecol->size(); ++
i){
203 int numGoodChannels = 0;
204 bool foundmaxcellx =
false;
205 bool foundmaxcelly =
false;
206 for(
unsigned int iPlane = 0; iPlane<
fgeom->
NPlanes();++iPlane){
210 foundmaxcellx =
true;
214 foundmaxcelly =
true;
216 for(
unsigned int iCell = 0;iCell<p->
Ncells();++iCell){
217 if(!(
fbadc->
IsBad(iPlane,iCell))){ ++numGoodChannels; }
228 for(
unsigned int iSlice = 0; iSlice < slicelist.
size(); iSlice++){
230 if(slicelist[iSlice]->IsNoise() || slicelist[iSlice]->NCell() >
fMaxHitCut){
continue; }
233 if(slicelist[iSlice]->NCell() < numGoodChannels*
fMaxSliceHits){
238 avgX = slicelist[iSlice]->MeanX();
239 avgY = slicelist[iSlice]->MeanY();
241 int minplanex = 10000;
245 std::vector<trk::LocatedChan> xcellchans;
246 xcellchans.
reserve(slicelist[iSlice]->NXCell());
247 double zxmin = 100000;
249 for(
unsigned int nx = 0; nx < slicelist[iSlice]->NXCell(); ++nx){
256 xchan.SetCenter(xyz[0],xyz[1],xyz[2]);
258 xcellchans.push_back(xchan);
259 if(xchan.Plane() < minplanex){ minplanex = xchan.Plane(); }
260 if(xchan.Plane() > maxplanex){ maxplanex = xchan.Plane(); }
265 if(xyz[2]>zxmax){ zxmax = xyz[2]; }
266 if(xyz[2]<zxmin){ zxmin = xyz[2]; }
272 unsigned int nxzbreak = zxmax-zxmin > 0? (
int)
floor((zxmax-zxmin)/100) : 0;
274 double dz = (zxmax-zxmin)/((
double)nxzbreak);
275 double xzbreaks[nxzbreak+2];
277 for(
unsigned int ibrk = 1; ibrk<=nxzbreak; ++ibrk){
278 xzbreaks[ibrk] = zxmin+((double)ibrk)*
dz;
280 xzbreaks[nxzbreak+1] = zxmax;
281 for(
unsigned int ibrk = 1; ibrk < nxzbreak+2; ++ibrk){
284 for(
unsigned int xchan = 0;xchan<xcellchans.size();++xchan){
286 xcellchans[xchan].GetCenter(xyz);
287 if(xyz[2] >= xzbreaks[ibrk-1] && xyz[2] <=xzbreaks[ibrk]){
292 std::vector<double> zrng;
293 zrng.push_back(xzbreaks[ibrk-1]);
294 zrng.push_back(xzbreaks[ibrk]);
296 if(ntot>0){
avgXPos.push_back(xtot/ntot); }
300 int minplaney = 10000;
304 std::vector<trk::LocatedChan> ycellchans;
305 double zymin = 100000;
307 for(
unsigned int ny = 0; ny < slicelist[iSlice]->NYCell(); ++ny){
314 ychan.SetCenter(xyz[0],xyz[1],xyz[2]);
316 ycellchans.push_back(ychan);
317 if(ychan.Plane() < minplaney){ minplaney = ychan.Plane(); }
318 if(ychan.Plane() > maxplaney){ maxplaney = ychan.Plane(); }
323 if(xyz[2]>zymax){ zymax = xyz[2]; }
324 if(xyz[2]<zymin){ zymin = xyz[2]; }
328 int nyzbreak = zymax-zymin > 0? (
int)
floor((zymax-zymin)/100) : 0;
330 dz = (zymax-zymin)/((
double)nyzbreak);
331 std::vector<double> yzbreaks;
332 yzbreaks.push_back(zymin);
333 for(
int ibrk = 1; ibrk<=nyzbreak; ++ibrk){
334 yzbreaks.push_back(zymin+((
double)ibrk)*dz);
336 yzbreaks.push_back(zymax);
337 for(
unsigned int ibrk = 1; ibrk < yzbreaks.size(); ++ibrk){
340 for(
unsigned int ychan = 0;ychan<ycellchans.size();++ychan){
342 ycellchans[ychan].GetCenter(xyz);
343 if(xyz[2] >= yzbreaks[ibrk-1] && xyz[2] <=yzbreaks[ibrk]){
348 std::vector<double> zrng;
349 zrng.push_back(yzbreaks[ibrk-1]);
350 zrng.push_back(yzbreaks[ibrk]);
352 if(ntot>0){
avgYPos.push_back(ytot/ntot); }
357 std::vector<rb::Track> xTracks;
367 std::vector<rb::Track> yTracks;
376 std::vector<rb::Track> alltracks;
377 for(
unsigned int i = 0;
i<xTracks.size();++
i){
378 xTracks[
i].SetID(iSlice);
379 alltracks.push_back(xTracks[
i]);
381 for(
unsigned int i = 0;
i<yTracks.size();++
i){
382 yTracks[
i].SetID(iSlice);
383 alltracks.push_back(yTracks[
i]);
385 for(std::vector<rb::Track>::const_iterator
it = alltracks.begin();
it != alltracks.end(); ++
it){
386 trackcol->push_back(*
it);
392 evt.
put(std::move(trackcol));
393 evt.
put(std::move(assns));
400 std::vector<rb::Track> outtracks;
417 std::vector<rb::Track> trackSegs;
418 std::vector< std::vector<double> > iP;
419 std::vector<bool> zProps;
420 std::vector< std::vector<bool> > trkHits;
421 std::vector< std::vector<std::array<double, 2> > > trkTrajs;
422 std::vector< std::vector<std::array<double, 2> > > trkDirs;
423 bool initTrack =
true;
426 int ntracks =
SegmentFinder(kgeo,iP,zProps,trkHits,trkTrajs,trkDirs);
429 std::vector< double > Chi2(trkHits.size(),0.0);
430 std::vector< int > nTrkHits(trkHits.size());
431 std::vector< bool > ignoreTrk(trkHits.size(),
false);
433 unsigned int propDirIters = 3;
435 for(
unsigned int iprop = 0; iprop < propDirIters; ++iprop){
438 int largestTrack = 0;
439 for(
unsigned int it = 0;
it< trkHits.size(); ++
it){
445 Chi2[
it] =
FilterTracker(kgeo,propDir,zProps[it],trkHits[it],trkTrajs[it],trkDirs[it],iP[it],nTrkHits[it]);
448 Chi2[
it] =
FilterTracker(kgeo,propDir,zProps[it],trkHits[it],trkTrajs[it],trkDirs[it],iP[it],nTrkHits[it]);
450 if(nTrkHits[it] > largestTrack){ largestTrack = nTrkHits[
it]; }
457 RemoveSameTracks(trkHits,trkTrajs,trkDirs,Chi2,iP,zProps,nTrkHits,ignoreTrk);
459 propDir = (-1*propDir);
467 for(
unsigned int i = 0;
i< trkHits.size();
i++){
470 if(nTrkHits[
i] >= minFiltHits && !ignoreTrk[
i]){
482 std::vector<bool> sortedhits = trkHits[
i];
486 trkHits[
i] = sortedhits;
495 int nhitsold = nTrkHits[
i];
501 while(filter && filterIter < maxFilter){
507 Chi2[
i] =
FilterOnly(trkHits[i],trkTrajs[i],trkDirs[i],iP[i],propDir,zProps[i],trkOutlier);
511 if(trkOutlier != -1){
512 trkHits[
i][trkOutlier] =
false;
519 if(nTrkHits[i] < minFiltHits){
525 if(nTrkHits[i] != nhitsold){
537 if(nhitsold == nTrkHits[i] && filterIter > 0){
541 if(iprop == propDirIters-1 || !initTrack){
544 std::vector<bool> newTrkHits;
545 CheckTrack(kgeo,trkHits[i],trkTrajs[i],zProps[i],newTrkHits);
551 if(nhitsnew > 0 && nTrkHits[i] > 0){
555 if(nTrkHits[i] < minFiltHits){ filter =
false; }
558 if(filter && zProps[i]){
568 if(nhitsnew >= minFiltHits){
570 zProps.push_back(zProps[i]);
572 trkHits.push_back(newTrkHits);
573 trkTrajs.push_back(trkTrajs[i]);
574 trkDirs.push_back(trkDirs[i]);
575 nTrkHits.push_back(nhitsnew);
576 ignoreTrk.push_back(
false);
595 nhitsold = nTrkHits[
i];
601 if(!initTrack){
break; }
605 unsigned int loopMinHits =
fMinHits;
607 for (
unsigned int s=0;
s<trkHits.size(); ++
s){
608 if ((.85*nTrkHits[
s]) > loopMinHits && !ignoreTrk[
s]) {
609 unsigned int trkMinHits = .85*nTrkHits[
s];
610 loopMinHits = trkMinHits;
617 if(trkHits.size() > 0){
620 double chiOverLength = 99999;
622 for(
unsigned int i = 0;
i<trkHits.size();
i++){
623 double nhits = nTrkHits[
i];
624 if(nhits < loopMinHits || ignoreTrk[
i]){
continue; }
625 double chiOverLengthnew = Chi2[
i]/nhits;
626 if(chiOverLengthnew < chiOverLength){
628 chiOverLength = chiOverLengthnew;
637 trkTrajs[bestTrack][firstHit][1],
638 trkDirs[bestTrack][firstHit][0],
639 trkDirs[bestTrack][firstHit][1]);
643 for(
int ihit = 0; ihit <
fNChans; ++ihit){
644 if(trkHits[bestTrack][ihit]){
645 if(zProps[bestTrack]){ tracknew.
Add(
fSliceChans[ihit].GetHit()); }
647 if(ihit != firstHit){ tracknew.
AppendTrajectoryPoint(trkTrajs[bestTrack][ihit][0],trkTrajs[bestTrack][ihit][1]); }
651 std::vector<int> extratrpts;
656 oldtrpt.Z() == tracknew.
TrajectoryPoint(itraj).Z()){ extratrpts.push_back(itraj);}
658 std::sort(extratrpts.begin(),extratrpts.end());
659 for(
int rmpt = ((
int)extratrpts.size()-1);rmpt>=0;--rmpt){
662 if(zProps[bestTrack]){
667 TVector3
dir = tracknew.
Dir();
668 double startslope = dir.X()/dir.Z();
669 if(
fView == 1){startslope = dir.Y()/dir.Z(); }
672 int startCell = 10000;
673 if(startslope < 0){ startCell = -1; }
674 for(
unsigned int iHit = 0; iHit<tracknew.
NCell(); ++iHit){
676 if(startslope<0 && tracknew.
Cell(iHit)->
Cell()>startCell){startCell = tracknew.
Cell(iHit)->
Cell();}
677 if(startslope>0 && tracknew.
Cell(iHit)->
Cell()<startCell){startCell = tracknew.
Cell(iHit)->
Cell();}
680 if(startCell != 10000){
681 std::array<double, 2> bestStart;
685 double deltaz = scell->
HalfD();
686 double deltav = scell->
HalfW();
688 if(
fView == 0){ kgeo.
BestTrackPoint(cpos[2],cpos[0],deltaz,deltav,bestStart,start.Z(),start.X(),dir.Z(),dir.X()); }
689 else if(
fView == 1){ kgeo.
BestTrackPoint(cpos[2],cpos[1],deltaz,deltav,bestStart,start.Z(),start.Y(),dir.Z(),dir.Y()); }
690 tracknew.
SetStart(bestStart[1],bestStart[0]);
694 TVector3 stop = tracknew.
Stop();
699 double stopslope = dir.X()/dir.Z();
700 if(
fView == 1){stopslope = dir.Y()/dir.Z(); }
703 int stopCell = 10000;
704 if(stopslope > 0){ stopCell = -1; }
705 for(
unsigned int iHit = 0; iHit<tracknew.
NCell(); ++iHit){
707 if(stopslope<0 && tracknew.
Cell(iHit)->
Cell()<stopCell){stopCell = tracknew.
Cell(iHit)->
Cell();}
708 if(stopslope>0 && tracknew.
Cell(iHit)->
Cell()>stopCell){stopCell = tracknew.
Cell(iHit)->
Cell();}
711 if(stopCell!= 10000){
712 std::array<double, 2> bestStop;
716 double deltaz = stcell->
HalfD();
717 double deltav = stcell->
HalfW();
719 if(
fView == 0){ kgeo.
BestTrackPoint(cpos[2],cpos[0],deltaz,deltav,bestStop,stop.Z(),stop.X(),dir.Z(),dir.X()); }
720 else if(
fView == 1){ kgeo.
BestTrackPoint(cpos[2],cpos[1],deltaz,deltav,bestStop,stop.Z(),stop.Y(),dir.Z(),dir.Y()); }
732 std::vector<TVector3> trjpoints;
737 reverse(trjpoints.begin(),trjpoints.end());
738 tracknew.
SetStart(trjpoints[0].
X(), trjpoints[0].
Z());
739 for(
unsigned int i = 1;
i<trjpoints.size(); ++
i){
747 std::vector<TVector3> trjpoints;
752 reverse(trjpoints.begin(),trjpoints.end());
753 tracknew.
SetStart(trjpoints[0].
Y(), trjpoints[0].
Z());
754 for(
unsigned int i = 1;
i<trjpoints.size(); ++
i){
763 outtracks.push_back(tracknew);
768 if(fNChans < (
int)
fMinHits){ Done =
true; }
773 if(nTimes == 2){ Done =
true; }
778 if(nTimes == 2){ Done =
true; }
788 std::vector< std::vector<double> > &iP,
789 std::vector< bool> &zProp,
790 std::vector< std::vector<bool> > &trkHits,
791 std::vector< std::vector<std::array<double, 2> > > &trkTrajs,
792 std::vector< std::vector<std::array<double, 2> > > &trkDirs)
796 std::vector<bool> hitsTemplate(
fNChans,
false);
797 std::vector<std::array<double, 2> > trajs(
fNChans);
799 std::vector<bool> usedVert(
fNChans,
false);
800 std::vector<bool> usedHor(
fNChans,
false);
805 int firsthitplane = 0;
806 SingleSegment(iP,zProp,trkHits,trkTrajs,trkDirs,firsthitplane);
815 firsthitplane = planeextent/2+
fSliceChans[0].Plane();
816 SingleSegment(iP,zProp,trkHits,trkTrajs,trkDirs,firsthitplane);
817 if((
int)iP.size() > ntracks){
820 if(ntracks < 2){ firstiHitPlane = firsthitplane; }
825 firsthitplane = firsthitplane + planeextent/4;
826 SingleSegment(iP,zProp,trkHits,trkTrajs,trkDirs,firsthitplane);
827 if((
int)iP.size() > ntracks){
830 if(ntracks < 2){ firstiHitPlane = firsthitplane; }
835 if(ntracks > 1){ firstiHitPlane = 0; }
843 for(
int iHit =
fNChans-1; iHit > 0; --iHit){
846 if(
fSliceChans[iHit].Plane() >= firstiHitPlane){
continue; }
849 for(
int fHit = iHit-1; fHit >= 0; --fHit){
856 int deltap = (iPlane-fPlane-2)/2;
859 if(deltap > maxgap){
continue; }
872 if(iPlane == fPlane && !usedVert[iHit]){
874 if(TMath::Abs(iCell-fCell)<=1){
877 zProp.push_back(
false);
879 std::vector<bool>
hits = hitsTemplate;
882 hits[celliHitid] =
true;
883 hits[cellfHitid] =
true;
884 trkHits.push_back(hits);
885 trajs[celliHitid][0] = ixyz[
fView];
886 trajs[celliHitid][1] = ixyz[2];
887 trajs[cellfHitid][0] = fxyz[
fView];
888 trajs[cellfHitid][1] = fxyz[2];
889 trkTrajs.push_back(trajs);
890 dirs[celliHitid][0] = 1.0;
891 dirs[celliHitid][1] = 0.0;
892 dirs[cellfHitid] = dirs[celliHitid];
893 trkDirs.push_back(dirs);
895 std::vector<double>
P(3);
900 usedVert[iHit] =
true;
904 else if(iCell == fCell && !usedHor[iHit]){
908 zProp.push_back(
true);
910 std::vector<double>
P(3);
913 P[1] = Pconst/deltaZ;
914 P[2] = 2*P[1]/deltaZ;
916 usedHor[iHit] =
true;
917 std::vector<bool>
hits = hitsTemplate;
920 trkHits.push_back(hits);
921 trajs[iHit][0] = ixyz[
fView];
922 trajs[iHit][1] = ixyz[2];
923 trajs[fHit][0] = fxyz[
fView];
924 trajs[fHit][1] = fxyz[2];
925 trkTrajs.push_back(trajs);
928 dirs[fHit] = dirs[iHit];
929 trkDirs.push_back(dirs);
932 else if(iPlane != fPlane && iCell != fCell){
934 double xyz[3] = {fxyz[0], fxyz[1], fxyz[2]+
fCellL-0.001};
935 double xyz1[3] = {ixyz[0], ixyz[1], ixyz[2]-
fCellL-0.001};
936 double deltaZ = ixyz[2]-fxyz[2]-2.0*
fCellL;
938 double m = (1.0/(deltaZ+4*
fCellL));
944 m*=TMath::Abs((delta-4*
fCellW));
946 else if(iCell<fCell){
950 m*=TMath::Abs((delta+4*
fCellW));
956 if(iPlane>(2+fPlane)){
959 if(numMissed <= maxgap){
960 double slope = delta/deltaZ;
961 double dxyz[3] = {0, 0, 1/
sqrt(1+slope*slope)};
964 std::vector<double>
P(3);
967 P[1] = Pconst/deltaZ;
968 P[2] = 2*P[1]/deltaZ;
969 std::vector<bool>
hits = hitsTemplate;
972 trajs[iHit][0] = ixyz[
fView];
973 trajs[iHit][1] = ixyz[2];
974 trajs[fHit][0] = fxyz[
fView];
975 trajs[fHit][1] = fxyz[2];
976 dirs[iHit][0] = dxyz[
fView];
977 dirs[iHit][1] = dxyz[2];
978 dirs[fHit] = dirs[iHit];
980 trkHits.push_back(hits);
981 zProp.push_back(
true);
983 trkTrajs.push_back(trajs);
984 trkDirs.push_back(dirs);
997 std::vector< bool> &zProp,
998 std::vector< std::vector<bool> > &trkHits,
999 std::vector< std::vector<std::array<double, 2> > > &trkTrajs,
1000 std::vector< std::vector<std::array<double, 2> > > &trkDirs,
1008 std::vector<double> vv;
1009 std::vector<double>
zz;
1010 std::vector<double>
w;
1014 double minz = 100000;
1016 double minv = 10000;
1017 double maxv = -10000;
1018 for(
int ihit = 0; ihit <
fNChans; ++ihit){
1019 if(channels[ihit].Plane() >= firsthitplane){
1022 channels[ihit].GetCenter(hitxyz);
1023 vv.push_back(hitxyz[
fView]);
1024 zz.push_back(hitxyz[2]);
1025 w.push_back(channels[ihit].GetHit()->PE());
1026 petot+=channels[ihit].GetHit()->PE();
1027 if(hitxyz[fView] > maxv){ maxv = hitxyz[
fView]; }
1028 else if (hitxyz[fView] < minv){ minv = hitxyz[
fView]; }
1029 if(hitxyz[2] > maxz){ maxz = hitxyz[2]; }
1030 else if(hitxyz[2] < minz){ minz = hitxyz[2]; }
1035 double linfitchi2 = 0.0;
1036 int nhits = vv.size();
1037 if(nhits < 3){
return 0.0; }
1038 for(
int i = 0;
i < nhits; ++
i){
1041 double z1[1],v1[1],z2[1],v2[1];
1044 double chindf = linfitchi2/(nhits-2);
1046 if(chindf > 3){
return chindf; }
1056 double deltav = v2[0] - v1[0];
1057 double deltaz = z2[0] - z1[0];
1058 double dirv = deltav/(deltav*deltav+deltaz*deltaz);
1059 double dirz = deltaz/(deltav*deltav+deltaz*deltaz);
1062 std::vector<std::array<double, 2> > trajs(fNChans);
1063 std::vector<std::array<double, 2> >
dirs(fNChans);
1065 for(
int ihit = 0; ihit <
fNChans; ++ihit){
1069 channels[ihit].GetCenter(hitxyz);
1071 trajs[ihit][1] = hitxyz[2];
1072 trajs[ihit][0] = (dirv/dirz)*(hitxyz[2]-z1[0])+v1[0];
1075 trajs[ihit][0] = hitxyz[
fView];
1076 trajs[ihit][1] = (dirz/dirv)*(hitxyz[
fView]-v1[0])+z1[0];
1078 dirs[ihit][0] = dirv;
1079 dirs[ihit][1] = dirz;
1084 double vdelta = maxv-minv;
1085 double zdelta = maxz-minz;
1087 double Pconst = vdelta*vdelta/12.0;
1088 std::vector<double>
P(3);
1090 P[1] = Pconst/zdelta;
1091 P[2] = 2*P[1]/zdelta;
1093 trkHits.push_back(hits);
1094 trkTrajs.push_back(trajs);
1095 trkDirs.push_back(dirs);
1096 zProp.push_back(prop);
1104 int zDir,
bool zProp, std::vector<bool> &trkHits,
1105 std::vector<std::array<double, 2> > &trkTrajs,
1106 std::vector<std::array<double, 2> > &trkDirs,
1107 std::vector<double> &errorCov,
int &outnhits,
1120 if(nHits == 0){
return 0.0; }
1125 if(zDir < 0){ iFirstHit =
fNChans-1; }
1127 bool firstHitfound =
false;
1128 while(!firstHitfound && iFirstHit<fNChans && iFirstHit >= 0){
1130 if(trkHits[iFirstHit]){ firstHitfound =
true; }
1131 else{iFirstHit+=zDir;}
1134 if(trkHits[iFirstHit]){ firstHitfound =
true; }
1135 else{ iFirstHit+=zDir; }
1139 else{ iFirstHit = firsthit; }
1141 unsigned short iPlane = channels[iFirstHit].Plane();
1142 unsigned short iCell = channels[iFirstHit].Cell();;
1145 int iCurrentHit = 0;
1146 if(zDir < 0){ iCurrentHit =
fNChans-1; }
1147 bool currentHitfound =
false;
1148 while(!currentHitfound && iCurrentHit<fNChans && iCurrentHit >= 0){
1150 if(channels[iCurrentHit].Plane() == iPlane){ currentHitfound =
true; }
1151 else{iCurrentHit+=zDir;}
1154 if(channels[iCurrentHit].Cell() == iCell){ currentHitfound =
true; }
1155 else{ iCurrentHit+=zDir; }
1160 const std::array<double, 2> &
start = trkTrajs[iFirstHit];
1161 std::array<double, 2>
dir = trkDirs[iFirstHit];
1162 int iCurrentTrackHit = iFirstHit;
1168 position = start[0];
1169 otherpos = start[1];
1170 if(dir[1] == 0){ dir[1] = 0.0001; }
1171 slope = dir[0]/dir[1];
1174 position = start[1];
1175 otherpos = start[0];
1176 if(dir[0] == 0){ dir[0] = 0.0001; }
1177 slope = dir[1]/dir[0];
1180 double R = cellVarW/3.0;
1181 if(!zProp){ R = cellVarL/3.0; };
1186 int iNextHit = iCurrentHit;
1195 else{ outnhits = nHits; }
1197 int nTrackHits = outnhits;
1207 short currentHitPlane = iPlane;
1208 int expectedHitPlane = currentHitPlane;
1209 short currentHitCell = iCell;
1210 short expectedHitCell = currentHitCell;
1212 std::array<double, 2> bestpos;
1214 while(((iCurrentHit < (
fNChans-1) && zDir > 0) || (iCurrentHit > 0 && zDir < 0)) && !done){
1216 short nextHitPlane = channels[iNextHit].Plane();
1217 short nextHitCell = channels[iNextHit].Cell();
1221 currentHitPlane = channels[iCurrentTrackHit].Plane();
1222 currentHitCell = channels[iCurrentTrackHit].Cell();
1226 if(zProp){sameLayerTest = ((nextHitPlane-currentHitPlane) == (expectedHitPlane-currentHitPlane)); }
1227 else{ sameLayerTest = (((nextHitCell-currentHitCell) == (expectedHitCell-currentHitCell))); }
1232 double peWeight = 0;
1233 double layerState[3] = {state[0], state[1], state[2]};
1234 double layerP[3] = {errorCov[0], errorCov[1], errorCov[2]};
1237 int firstLayerHit = iNextHit;
1239 if(zProp){
while(k < fNChans && k >= 0 && channels[k].Plane() == nextHitPlane){ k+=zDir; } }
1240 else{
while(k < fNChans && k >= 0 && channels[k].Cell() == nextHitCell){ k+=zDir; } }
1241 int nadd = zDir*(k-iNextHit);
1242 bool layerTrackHits[nadd];
1243 int iLayerHit = firstLayerHit;
1244 int fLayerHit = k-zDir;
1245 for(
int i = 0;
i < nadd; ++
i){
1246 layerTrackHits[
i] =
false;
1247 if(trkHits[iNextHit]){
1248 trkHits[iNextHit] =
false;
1251 iCurrentHit = iNextHit;
1252 iNextHit = iCurrentHit + zDir;
1254 bool doneAdding =
false;
1255 if(nadd == 0){ doneAdding =
true;}
1258 int layerWindowlow = 0;
1259 int layerWindowhigh = 999999;
1261 if(((state[2] <= 0 && zDir < 0) || (state[2] >= 0 && zDir > 0))){
1263 layerWindowlow = channels[iCurrentTrackHit].Plane();
1264 if(
abs(nextHitCell-currentHitCell)>1){
1265 layerWindowlow = channels[iCurrentTrackHit].Plane()-2;
1266 layerWindowhigh = channels[iCurrentTrackHit].Plane()+6;
1272 layerWindowhigh = channels[iCurrentTrackHit].Plane();
1273 if(
abs(nextHitCell-currentHitCell)>1){
1274 layerWindowhigh = channels[iCurrentTrackHit].Plane()+2;
1275 layerWindowlow = channels[iCurrentTrackHit].Plane()-6;
1281 double minDistBest = -1;
1282 int bestHitIndex = 9999;
1283 if(nadd == 1 && !layerTrackHits[0]){ bestHitIndex = 0; }
1285 for(
int i = 0;
i<nadd;++
i){
1287 if(!layerTrackHits[
i]){
1288 int currChanIdx = firstLayerHit+i*zDir;
1291 if(zProp && (channels[currChanIdx].Cell()<layerWindowlow || channels[currChanIdx].Cell()>layerWindowhigh)){
continue; }
1292 else if(!zProp && (channels[currChanIdx].Plane()<layerWindowlow || channels[currChanIdx].Plane()>layerWindowhigh)){
continue; }
1296 channels[currChanIdx].GetCenter(xyzHit);
1297 if(zProp){kgeo.
BestTrackPoint(xyzHit[2],xyzHit[
fView],fCellL,fCellW,bestpos,layerState[1],layerState[0],dir[1],dir[0]);}
1298 else{kgeo.
BestTrackPoint(xyzHit[2],xyzHit[
fView],fCellL,fCellW,bestpos,layerState[0],layerState[1],dir[1],dir[0]);}
1300 double minHitDist = 0.0;
1301 if(zProp){ minHitDist =
abs(bestpos[1]-(layerState[0]+layerState[2]*(bestpos[0]-layerState[1])));}
1302 else{ minHitDist =
abs(bestpos[0]-(layerState[0]+layerState[2]*(bestpos[1]-layerState[1]))); }
1303 if((minHitDist < minDistBest) || minDistBest < 0){
1304 bestHitIndex = i*zDir;
1305 minDistBest = minHitDist;
1310 if(bestHitIndex != 9999){
1311 int bestChanIdx = firstLayerHit+bestHitIndex;
1313 channels[bestChanIdx].GetCenter(xyzHit);
1316 if(zProp){ delta = xyzHit[2] - layerState[1]; }
1317 else{ delta = xyzHit[
fView] - layerState[1]; }
1318 Q[0] = delta*delta*Q[2] + layerP[0];
1322 xMinus[0] = layerState[0] + delta*layerState[2];
1323 xMinus[2] = layerState[2];
1324 pMinus[0] = layerP[0] + 2.0*delta*layerP[1] + delta*delta*layerP[2] + Q[0];
1325 pMinus[1] = layerP[1] + delta*layerP[2] + Q[1];
1326 pMinus[2] = layerP[2] + Q[2];
1329 double denom = (pMinus[0]+
R);
1330 K[0] = pMinus[0]/denom;
1331 K[1] = pMinus[1]/denom;
1334 double candidatePos;
1337 candidatePos = xyzHit[
fView];
1342 candidatePos = xyzHit[2];
1343 ipos = xyzHit[
fView];
1348 xPredict[0] = xMinus[0]+K[0]*(candidatePos-xMinus[0]);
1349 xPredict[2] = xMinus[2]+K[1]*(candidatePos-xMinus[0]);
1350 cPredict[0] = pMinus[0]*R/denom;
1351 cPredict[1] = pMinus[1]*R/denom;
1352 cPredict[2] = (pMinus[2]*pMinus[0]-pMinus[1]*pMinus[1]+pMinus[2]*
R)/denom;
1354 double r2 = (candidatePos - xPredict[0])*(candidatePos - xPredict[0]);
1355 double sl2 = xPredict[2]*xPredict[2];
1356 double deltaChi = 0.0;
1358 deltaChi = (r2*(1+sl2)*(1+sl2))/
1359 (((1+sl2)*(1+sl2)*(cellVarW/3+sl2*cellVarL/3+cPredict[0]))+r2*sl2*cPredict[2]/2);
1362 deltaChi = (r2*(1+sl2)*(1+sl2))/
1363 (((1+sl2)*(1+sl2)*(cellVarL/3+sl2*cellVarW/3+cPredict[0]))+r2*sl2*cPredict[2]/2);
1366 if(deltaChi <= maxDeltaChi){
1369 trkHits[bestChanIdx] =
true;
1372 layerTrackHits[(bestHitIndex/zDir)] =
true;
1373 currentHitPlane = channels[bestChanIdx].Plane();
1374 currentHitCell = channels[bestChanIdx].Cell();
1375 float currentPE = channels[bestChanIdx].GetHit()->PE();
1376 iCurrentTrackHit = bestChanIdx;
1378 layerState[0] = xPredict[0];
1379 layerState[1] = xPredict[1];
1380 layerState[2] = xPredict[2];
1381 layerP[0] = cPredict[0];
1382 layerP[1] = cPredict[1];
1383 layerP[2] = cPredict[2];
1384 peWeight = currentPE;
1387 layerState[0] = (layerState[0]*peWeight+xPredict[0]*currentPE)/(peWeight+currentPE);
1388 layerState[1] = (layerState[1]*peWeight+xPredict[1]*currentPE)/(peWeight+currentPE);
1389 layerState[2] = (layerState[2]*peWeight+xPredict[2]*currentPE)/(peWeight+currentPE);
1390 layerP[0] = (layerP[0]*peWeight+cPredict[0]*currentPE)/(peWeight+currentPE);
1391 layerP[1] = (layerP[1]*peWeight+cPredict[1]*currentPE)/(peWeight+currentPE);
1392 layerP[2] = (layerP[2]*peWeight+cPredict[2]*currentPE)/(peWeight+currentPE);
1393 peWeight+=currentPE;
1397 dir[0] = layerState[2]/
sqrt(1.0+layerState[2]*layerState[2]);
1398 dir[1] = 1.0/
sqrt(1.0+layerState[2]*layerState[2]);
1400 int newLayerWindowlow = currentHitCell-
fGapAllowed;
1401 int newLayerWindowhigh = currentHitCell+
fGapAllowed;
1403 while(newLayerWindowlow > 0 &&
1404 (
fbadc->
IsBad(currentHitPlane,newLayerWindowlow)
1406 --newLayerWindowlow;
1409 (
fbadc->
IsBad(currentHitPlane,newLayerWindowhigh)
1411 ++newLayerWindowhigh;
1414 layerWindowlow = newLayerWindowlow;
1415 layerWindowhigh = newLayerWindowhigh;
1418 if((currentHitCell - layerWindowlow) < (layerWindowhigh - currentHitCell)){
1422 int newLayerWindowlow = currentHitCell-
fGapAllowed;
1423 while(newLayerWindowlow > 0 &&
1424 (
fbadc->
IsBad(currentHitPlane,newLayerWindowlow)
1426 --newLayerWindowlow;
1428 layerWindowlow = newLayerWindowlow;
1435 int newLayerWindowhigh = currentHitCell+
fGapAllowed;
1437 (
fbadc->
IsBad(currentHitPlane,newLayerWindowhigh)
1439 ++newLayerWindowhigh;
1441 layerWindowhigh = newLayerWindowhigh;
1448 dir[0] = 1.0/
sqrt(1.0+layerState[2]*layerState[2]);
1449 dir[1] = layerState[2]/
sqrt(1.0+layerState[2]*layerState[2]);
1451 if(((state[2] <= 0 && zDir < 0) || (state[2] >= 0 && zDir > 0))){
1452 int newLayerWindowlow = currentHitPlane;
1454 layerWindowlow = newLayerWindowlow;
1458 int newlayerWindowhightest = newLayerWindowhigh;
1460 (
fbadc->
IsBad(newlayerWindowhightest,currentHitCell)
1467 layerWindowhigh = newLayerWindowhigh;
1471 int newLayerWindowhigh = currentHitPlane;
1472 layerWindowhigh = newLayerWindowhigh;
1476 int newlayerWindowlowtest = newLayerWindowlow;
1478 (
fbadc->
IsBad(newlayerWindowlowtest,currentHitCell)
1485 layerWindowlow = newLayerWindowlow;
1489 if(((state[2] <= 0 && zDir < 0) || (state[2] >= 0 && zDir > 0))){
1490 layerWindowlow = currentHitPlane;
1499 int newLayerWindowhightest = newLayerWindowhigh;
1501 (
fbadc->
IsBad(newLayerWindowhightest,currentHitCell)
1507 layerWindowhigh = newLayerWindowhigh;
1511 layerWindowhigh = currentHitPlane;
1520 int newLayerWindowlowtest = newLayerWindowlow;
1522 (
fbadc->
IsBad(newLayerWindowlowtest,currentHitCell)
1528 layerWindowlow = newLayerWindowlow;
1534 else{ layerTrackHits[(bestHitIndex/zDir)] =
true; }
1537 else{ doneAdding =
true;}
1543 state[0] = layerState[0];
1544 state[1] = layerState[1];
1545 state[2] = layerState[2];
1546 errorCov[0] = layerP[0];
1547 errorCov[1] = layerP[1];
1548 errorCov[2] = layerP[2];
1552 if(fLayerHit < iLayerHit){
std::swap(iLayerHit,fLayerHit); }
1553 for(
int ihit = iLayerHit; ihit<=fLayerHit; ++ihit){
1555 trkTrajs[ihit][0] = state[0];
1556 trkTrajs[ihit][1] = state[1];
1557 trkDirs[ihit][0] = dir[0];
1558 trkDirs[ihit][1] = dir[1];
1565 if(fLayerHit < iLayerHit){
std::swap(iLayerHit,fLayerHit); }
1566 for(
int ihit = iLayerHit; ihit<=fLayerHit; ++ihit){
1568 trkTrajs[ihit][0] = state[1];
1569 trkTrajs[ihit][1] = state[0];
1570 trkDirs[ihit][0] = dir[0];
1571 trkDirs[ihit][1] = dir[1];
1574 expectedHitCell = currentHitCell+zDir;
1575 if(expectedHitCell < 0 || expectedHitCell >
NCells[
fView] ){ done =
true; }
1577 if(iNextHit > (
fNChans-1) || iNextHit < 0){ done =
true; }
1583 if(iNextHit > (
fNChans-1) || iNextHit < 0){ done =
true; }
1585 short nextHitPlane = channels[iNextHit].Plane();
1586 short nextHitCell = channels[iNextHit].Cell();
1591 channels[iNextHit].GetCenter(nextpt);
1594 double offset = zDir*(fCellL+2.0);
1595 currtraj[1] = state[1] +
offset;
1596 currtraj[0] = state[2]*offset + state[0];
1599 nextpt[
fView] = state[2]*(nextpt[2]-currtraj[1])+currtraj[0];
1606 double offset = zDir*(fCellW+2.0);
1607 currtraj[0] = state[1] +
offset;
1608 currtraj[1] = state[2]*offset + state[0];
1611 nextpt[2] = state[2]*(nextpt[
fView]-currtraj[0])+currtraj[1];
1612 if(nextpt[2] < 0 || nextpt[2] >
fgeom->
DetLength()){ done =
true; }
1615 if(misses > maxMisses){done =
true;}
1617 expectedHitCell = nextHitCell;
1618 expectedHitPlane = nextHitPlane;
1628 channels[iNextHit].GetCenter(nextpt);
1631 double offset = zDir*(fCellL+2.0);
1632 currtraj[1] = state[1] +
offset;
1633 currtraj[0] = state[2]*offset + state[0];
1636 nextpt[
fView] = state[2]*(nextpt[2]-currtraj[1])+currtraj[0];
1643 double offset = zDir*(fCellW+2.0);
1644 currtraj[0] = state[1] +
offset;
1645 currtraj[1] = state[2]*offset + state[0];
1648 nextpt[2] = state[2]*(nextpt[
fView]-currtraj[0])+currtraj[1];
1649 if(nextpt[2] < 0 || nextpt[2] >
fgeom->
DetLength()){ done =
true; }
1652 if(misses > maxMisses){done =
true;}
1654 expectedHitPlane = nextHitPlane;
1655 expectedHitCell = nextHitCell;
1667 std::vector< std::vector<std::array<double, 2> > > &trkTrajs,
1668 std::vector< std::vector<std::array<double, 2> > >&trkDirs,std::vector<double> &
chi2,
1669 std::vector< std::vector<double> > &pfiltState, std::vector<bool> &zProps,
1670 std::vector<int> &nTrkHits, std::vector<bool> &ignoreTrk)
1672 int nTracks = trkHits.size();
1673 if(nTracks == 0)
return;
1675 bool sameTrack[nTracks];
1676 for(
int st = 0; st < nTracks; ++st) sameTrack[st] =
false;
1679 for(
int iTrack = 0; iTrack < nTracks; iTrack++){
1681 bool inTrackDone =
false;
1682 if(sameTrack[iTrack] || ignoreTrk[iTrack]){inTrackDone =
true; }
1684 int cTrack = iTrack+1;
1686 while(cTrack < nTracks && !inTrackDone){
1687 if(!sameTrack[cTrack] && !ignoreTrk[cTrack]){
1691 if(zProps[iTrack] == zProps[cTrack] && nTrkHits[iTrack] == nTrkHits[cTrack]){
1693 for(
int ihit = 0; ihit <
fNChans; ++ihit){
1694 if(trkHits[iTrack][ihit] != trkHits[cTrack][ihit]){
1700 else{ same =
false; }
1703 if(chi2[iTrack] > chi2[cTrack]){
1704 sameTrack[iTrack] =
true;
1705 ignoreTrk[iTrack] =
true;
1710 sameTrack[cTrack] =
true;
1711 ignoreTrk[cTrack] =
true;
1719 if(sameTrack[iTrack]){ ++nSame; }
1730 std::vector<trk::LocatedChan> newfSignalChan;
1736 bool inTrack =
false;
1737 for(
unsigned int j = 0;
j < track.
NCell(); ++
j){
1740 if(chanhit == thit){
1745 if(!inTrack){ newfSignalChan.push_back(
fSliceChans[
i]); }
1747 fNChans = newfSignalChan.size();
1755 std::vector<std::array<double, 2> > &trkDirs, std::vector<double> &pEstimate,
1756 int zDir,
bool zProp,
int &newOutlierPos)
1769 const std::array<double, 2> &
start = trkTrajs[firstHit];
1770 std::array<double, 2>
dir = trkDirs[firstHit];
1773 double trajPointsFor[nhits][3];
1774 double pFor[nhits][4];
1775 double trajPointsBack[nhits][3];
1776 double pBack[nhits][4];
1783 otherpos = start[1];
1784 if(dir[1] == 0){ dir[1] = 0.0001; }
1785 slope = dir[0]/dir[1];
1789 otherpos = start[0];
1790 if(dir[0] == 0){ dir[0] = 0.0001; }
1791 slope = dir[1]/dir[0];
1794 std::vector<const rb::CellHit * > trackHitsNew(nhits);
1795 double hitCentersNew[nhits][3];
1799 trackHitsNew[ihit] = channels[
i].GetHit().get();
1800 channels[
i].GetCenter(hitCentersNew[ihit]);
1806 if(!zProp && (trackHitsNew[0]->Plane() == trackHitsNew[nhits-1]->Plane())){
1811 if(zDir > 0){ trkDirs[
i][0] = 1.0; }
1812 else{ trkDirs[
i][0] = -1.0; }
1813 trkDirs[
i][1] = 0.0;
1814 trkTrajs[
i][0] = hitCentersNew[ihit][
fView];
1815 trkTrajs[
i][1] = hitCentersNew[ihit][2];
1823 double R = cellVarW/3.0;
1824 if(!zProp){ R = cellVarL/3.0; }
1829 int iCurrentHit = 0;
1834 while(iCurrentHit<nhits && iNextHit<nhits){
1835 int searchPlane = trackHitsNew[iCurrentHit]->Plane();
1836 int searchCell = trackHitsNew[iCurrentHit]->Cell();
1840 for(
int i = iCurrentHit;
i<nhits;++
i){
1841 if(trackHitsNew[
i]->Plane() == searchPlane){
1844 peTot+= trackHitsNew[
i]->PE();
1850 for(
int i = iCurrentHit;
i<nhits; ++
i){
1851 if(trackHitsNew[
i]->Cell() == searchCell){
1854 peTot+= trackHitsNew[
i]->PE();
1860 double otherposOld = otherpos;
1861 double posOld =
pos;
1862 double slopeOld = slope;
1866 const double pEstimateOld[3] = { pEstimate[0], pEstimate[1], pEstimate[2] };
1871 while(iCurrentHit<iNextHit){
1874 z[0] = hitCentersNew[iCurrentHit][
fView];
1875 z[1] = hitCentersNew[iCurrentHit][2];
1878 z[0] = hitCentersNew[iCurrentHit][2];
1879 z[1] = hitCentersNew[iCurrentHit][
fView];
1881 delta = z[1] - otherposOld;
1884 Q[0] = delta*delta*Q[2];
1888 xMinus[0] = posOld + delta*slopeOld;
1889 xMinus[2] = slopeOld;
1890 pMinus[0] = pEstimateOld[0] + 2.0*delta*pEstimateOld[1] + delta*delta*pEstimateOld[2] + Q[0];
1891 pMinus[1] = pEstimateOld[1] + delta*pEstimateOld[2] + Q[1];
1892 pMinus[2] = pEstimateOld[2] + Q[2];
1895 double denom = (pMinus[0]+
R);
1896 K[0] = pMinus[0]/denom;
1897 K[1] = pMinus[1]/denom;
1899 const double weight = trackHitsNew[iCurrentHit]->PE()/peTot;
1900 pos+= weight*(xMinus[0] + K[0]*(z[0]-xMinus[0]));
1901 otherpos+= weight*z[1];
1903 slope+= weight*(xMinus[2] + K[1]*(z[0]-xMinus[0]));
1904 pEstimate[0]+= weight*(pMinus[0]*R/denom);
1905 pEstimate[1]+= weight*(pMinus[1]*R/denom);
1906 pEstimate[2]+= weight*((pMinus[2]*pMinus[0]-pMinus[1]*pMinus[1]+pMinus[2]*
R)/denom);
1909 while(iTraj<iNextHit){
1912 trajPointsFor[iTraj][0] =
pos;
1913 trajPointsFor[iTraj][1] = otherpos;
1916 trajPointsFor[iTraj][0] = otherpos;
1917 trajPointsFor[iTraj][1] =
pos;
1919 trajPointsFor[iTraj][2] = slope;
1920 pFor[iTraj][0] = pEstimate[0];
1921 pFor[iTraj][3] = pEstimate[2];
1924 iCurrentHit = iNextHit;
1929 iCurrentHit = nhits-1;
1933 while(iCurrentHit>=0 && iNextHit>=0){
1934 int searchPlane = trackHitsNew[iCurrentHit]->Plane();
1935 int searchCell = trackHitsNew[iCurrentHit]->Cell();
1939 for(
int i = iCurrentHit;
i>=0;--
i){
1940 if(trackHitsNew[
i]->Plane() == searchPlane){
1943 peTot+= trackHitsNew[
i]->PE();
1949 for(
int i = iCurrentHit;
i>=0; --
i){
1950 if(trackHitsNew[
i]->Cell() == searchCell){
1953 peTot+= trackHitsNew[
i]->PE();
1959 double otherposOld = otherpos;
1960 double posOld =
pos;
1961 double slopeOld = slope;
1965 const double pEstimateOld[3] = { pEstimate[0], pEstimate[1], pEstimate[2] };
1970 double weight = 1.0/((double)nhitsAdded);
1971 while(iCurrentHit>iNextHit){
1974 z[0] = hitCentersNew[iCurrentHit][
fView];
1975 z[1] = hitCentersNew[iCurrentHit][2];
1978 z[0] = hitCentersNew[iCurrentHit][2];
1979 z[1] = hitCentersNew[iCurrentHit][
fView];
1981 delta = z[1] - otherposOld;
1984 Q[0] = delta*delta*Q[2];
1988 xMinus[0] = posOld + delta*slopeOld;
1989 xMinus[2] = slopeOld;
1990 pMinus[0] = pEstimateOld[0] + 2.0*delta*pEstimateOld[1] + delta*delta*pEstimateOld[2] + Q[0];
1991 pMinus[1] = pEstimateOld[1] + delta*pEstimateOld[2] + Q[1];
1992 pMinus[2] = pEstimateOld[2] + Q[2];
1995 double denom = (pMinus[0]+
R);
1996 K[0] = pMinus[0]/denom;
1997 K[1] = pMinus[1]/denom;
1999 weight = trackHitsNew[iCurrentHit]->PE()/peTot;
2000 pos+= weight*(xMinus[0] + K[0]*(z[0]-xMinus[0]));
2001 otherpos+= weight*z[1];
2003 slope+= weight*(xMinus[2] + K[1]*(z[0]-xMinus[0]));
2004 pEstimate[0]+= weight*(pMinus[0]*R/denom);
2005 pEstimate[1]+= weight*(pMinus[1]*R/denom);
2006 pEstimate[2]+= weight*((pMinus[2]*pMinus[0]-pMinus[1]*pMinus[1]+pMinus[2]*
R)/denom);
2010 while(iTraj>iNextHit){
2012 trajPointsBack[iTraj][0] =
pos;
2013 trajPointsBack[iTraj][1] = otherpos;
2016 trajPointsBack[iTraj][0] = otherpos;
2017 trajPointsBack[iTraj][1] =
pos;
2019 trajPointsBack[iTraj][2] = slope;
2020 pBack[iTraj][0] = pEstimate[0];
2021 pBack[iTraj][3] = pEstimate[2];
2024 iCurrentHit = iNextHit;
2029 double worstchi = 0;
2033 double trajWeight[2] = {0.5, 0.5};
2034 double dirWeight[2] = {0.5, 0.5};
2037 double deltaChi = 0.0;
2039 if(trackHitsNew[ihit]->Plane() == trackHitsNew[0]->Plane()){
2040 trajWeight[0] = 0.0, trajWeight[1] = 1.0;
2041 dirWeight[0] = 0.0, dirWeight[1] = 1.0;
2042 p0 = pBack[ihit][0];
2043 p3 = pBack[ihit][3];
2045 else if(trackHitsNew[ihit]->Plane() == trackHitsNew[(nhits-1)]->Plane()){
2046 trajWeight[0] = 1.0, trajWeight[1] = 0.0;
2047 dirWeight[0] = 1.0, dirWeight[1] = 0.0;
2052 double trajdenom = 1.0/
sqrt(pFor[ihit][0]) + 1.0/
sqrt(pBack[ihit][0]);
2053 trajWeight[0] = (1.0/
sqrt(pFor[ihit][0]))/trajdenom, trajWeight[1] = (1.0/
sqrt(pBack[ihit][0]))/trajdenom;
2054 double dirdenom = 1.0/
sqrt(pFor[ihit][3]) + 1.0/
sqrt(pBack[ihit][3]);
2055 dirWeight[0] = (1.0/
sqrt(pFor[ihit][3]))/dirdenom, dirWeight[1] = (1.0/
sqrt(pBack[ihit][3]))/dirdenom;
2056 p0 = 2.0/(trajdenom*trajdenom);
2057 p3 = 2.0/(dirdenom*dirdenom);
2059 trkTrajs[
i][0] = trajWeight[0]*trajPointsFor[ihit][0] + trajWeight[1]*trajPointsBack[ihit][0];
2060 trkTrajs[
i][1] = 0.5*trajPointsFor[ihit][1] + 0.5*trajPointsBack[ihit][1];
2061 slope = dirWeight[0]*trajPointsFor[ihit][2] + dirWeight[1]*trajPointsBack[ihit][2];
2062 double sl2 = slope*slope;
2063 trkDirs[
i][0] = slope/
sqrt(1.0+sl2);
2064 trkDirs[
i][1] = 1.0/
sqrt(1.0+sl2);
2065 double r2 = (hitCentersNew[ihit][
fView]-trkTrajs[
i][0])*(hitCentersNew[ihit][
fView]-trkTrajs[i][0]);
2066 deltaChi+= (r2*(1.0+sl2)*(1.0+sl2))/(((1.0+sl2)*(1.0+sl2)*(cellVarW/3.0+sl2*cellVarL/3.0+p0))+r2*sl2*p3/2.0);
2069 if(trackHitsNew[ihit]->Cell() == trackHitsNew[0]->Cell()){
2070 trajWeight[0] = 0.0, trajWeight[1] = 1.0;
2071 dirWeight[0] = 0.0, dirWeight[1] = 1.0;
2072 p0 = pBack[ihit][0];
2073 p3 = pBack[ihit][3];
2075 else if(trackHitsNew[ihit]->Cell() == trackHitsNew[(nhits-1)]->Cell()){
2076 trajWeight[0] = 1.0, trajWeight[1] = 0.0;
2077 dirWeight[0] = 1.0, dirWeight[1] = 0.0;
2082 double trajdenom = 1.0/
sqrt(pFor[ihit][0]) + 1.0/
sqrt(pBack[ihit][0]);
2083 trajWeight[0] = (1.0/
sqrt(pFor[ihit][0]))/trajdenom, trajWeight[1] = (1.0/
sqrt(pBack[ihit][0]))/trajdenom;
2084 double dirdenom = 1.0/
sqrt(pFor[ihit][3]) + 1.0/
sqrt(pBack[ihit][3]);
2085 dirWeight[0] = (1.0/
sqrt(pFor[ihit][3]))/dirdenom, dirWeight[1] = (1.0/
sqrt(pBack[ihit][3]))/dirdenom;
2086 p0 = 2.0/(trajdenom*trajdenom);
2087 p3 = 2.0/(dirdenom*dirdenom);
2089 trkTrajs[
i][0] = 0.5*trajPointsFor[ihit][0] + 0.5*trajPointsBack[ihit][0];
2090 trkTrajs[
i][1] = trajWeight[0]*trajPointsFor[ihit][1] + trajWeight[1]*trajPointsBack[ihit][1];
2091 slope = dirWeight[0]*trajPointsFor[ihit][2] + dirWeight[1]*trajPointsBack[ihit][2];
2092 double sl2 = slope*slope;
2093 trkDirs[
i][0] = 1.0/
sqrt(1.0+sl2);
2094 trkDirs[
i][1] = slope/
sqrt(1.0+sl2);
2095 double r2 = (hitCentersNew[ihit][2]-trkTrajs[
i][1])*(hitCentersNew[ihit][2]-trkTrajs[i][1]);
2096 deltaChi+= (r2*(1+sl2)*(1+sl2))/(((1+sl2)*(1+sl2)*(cellVarL/3.0+sl2*cellVarW/3.0+p0))+r2*sl2*p3/2.0);
2101 if(deltaChi>worstchi){
2102 worstchi = deltaChi;
2116 std::vector<bool> &trkHits, std::vector<std::array<double, 2> > trkTrajs,
2117 bool zProp, std::vector<bool> &newTrkHits)
2128 std::vector<std::array<double, 2> > trjPointsNew(nhits);
2129 std::vector<int> trkHitsIdx(nhits);
2133 trjPointsNew[ihit] = trkTrajs[
i];
2134 trkHitsIdx[ihit] =
i;
2139 unsigned int lastPlane = 10000;
2140 unsigned int lastCell = 10000;
2141 unsigned int lastHitPlane = 10000;
2142 unsigned int lastHitCell = 10000;
2143 bool posSlope =
true;
2144 bool spillover =
false;
2145 int spillovergap = 0;
2146 double spilloverprob = 1.0;
2151 for(
int itraj = 0; itraj < (nhits-1); ++itraj){
2154 if(trjPointsNew[itraj][0] == trjPointsNew[itraj+1][0] && trjPointsNew[itraj][1] == trjPointsNew[itraj+1][1]){
continue; }
2157 std::vector<geo::OfflineChan> missedCells;
2162 trjPointsNew[itraj+1][0],trjPointsNew[itraj+1][1],missedCells,distance);
2164 if(totMissed == 0){
continue; }
2167 if(trjPointsNew[itraj][1] != trjPointsNew[itraj+1][1] &&
2168 ((trjPointsNew[itraj+1][0]-trjPointsNew[itraj][0])/(trjPointsNew[itraj+1][1]-trjPointsNew[itraj][1])) < 0){ sort =
true; }
2173 std::vector<geo::OfflineChan> orderedMissedCells(totMissed);
2174 std::vector<double> orderedDistance(totMissed);
2175 unsigned int firstCellIndex = 0;
2176 unsigned int lastCellIndex = 0;
2177 unsigned int fillCellIndex = 0;
2178 unsigned int plane = missedCells[0].Plane();
2179 for(
int iCell = 0; iCell<totMissed;++iCell){
2180 if(missedCells[iCell].Plane() != (
int)plane || (
int)iCell == totMissed-1){
2181 if(iCell == totMissed-1 && missedCells[iCell].Plane() == (
int)plane){lastCellIndex = iCell;}
2183 for(
int fillCell = (
int) lastCellIndex; fillCell>= (
int) firstCellIndex;--fillCell){
2184 orderedMissedCells[fillCellIndex] = missedCells[fillCell];
2185 orderedDistance[fillCellIndex] = distance[fillCell];
2188 if(iCell == totMissed-1 && missedCells[iCell].Plane() != (
int)plane){
2189 orderedMissedCells[totMissed-1] = missedCells[totMissed-1];
2190 orderedDistance[totMissed-1] = distance[totMissed-1];
2191 lastCellIndex = iCell;
2192 if(missedCells[iCell].Plane() !=
plane){firstCellIndex = iCell;}
2195 firstCellIndex = iCell;
2196 lastCellIndex = iCell;
2197 plane = missedCells[iCell].Plane();
2199 else{ lastCellIndex = iCell; }
2201 missedCells = orderedMissedCells;
2202 distance = orderedDistance;
2206 double currprob = 1.0;
2210 currentgap+=spillovergap;
2211 currprob = spilloverprob;
2215 for(
int iCell = 0; iCell < totMissed; ++iCell){
2217 bool ontrack =
false;
2218 bool inbadc =
false;
2221 for(
int ihit = 0; ihit < nhits; ++ihit){
2222 if(channels[trkHitsIdx[ihit]].Plane() == missedCells[iCell].Plane() &&
2223 channels[trkHitsIdx[ihit]].Cell() == missedCells[iCell].Cell()){
2225 lastHitPlane = missedCells[iCell].Plane();
2226 lastHitCell = missedCells[iCell].Cell();
2233 if(
fbadc->
IsBad(missedCells[iCell].Plane(),missedCells[iCell].Cell()) ||
2234 fbadc->
IsLowEfficiency(missedCells[iCell].Plane(),missedCells[iCell].Cell()) ){ inbadc =
true; }
2242 spilloverprob = 1.0;
2245 else if(!ontrack && !inbadc){
2247 double lenInCell = distance[iCell];
2248 double zave = 0.5*trjPointsNew[itraj][1]+0.5*trjPointsNew[itraj+1][1];
2254 bool foundrng =
false;
2275 if(!foundrng){ w =
avgYPos[0]; }
2279 bool foundrng =
false;
2300 if(!foundrng){ w =
avgXPos[0]; }
2313 if(missedCells[iCell].Plane() != lastPlane ||
2314 missedCells[iCell].Cell() != lastCell){
2323 if(iCell == (totMissed-1)){
2324 spillovergap = currentgap;
2325 spilloverprob = currprob;
2326 lastPlane = missedCells[iCell].Plane();
2327 lastCell = missedCells[iCell].Cell();
2334 if(currprob < minProb && currentgap > 1){
2335 if(lastHitPlane == 10000){
2336 lastHitPlane = missedCells[iCell].Plane();
2337 lastHitCell = missedCells[iCell].Cell();
2344 if(currentgap > biggestGap){ biggestGap = currentgap; }
2346 if(currprob < minProb && currentgap > 1){
2349 newTrkHits = trkHits;
2351 for(
int ihit = 0; ihit <
fNChans; ++ihit){
2353 if(!trkHits[ihit]){
continue; }
2355 if(channels[ihit].Plane() > lastHitPlane){
2356 trkHits[ihit] =
false;
2360 if(channels[ihit].Plane() < lastHitPlane){ newTrkHits[ihit] =
false; }
2362 if(channels[ihit].Plane() == lastHitPlane){
2363 if(posSlope && channels[ihit].Cell() > lastHitCell){
2364 trkHits[ihit] =
false;
2367 else if(posSlope && channels[ihit].Cell() <= lastHitCell){ newTrkHits[ihit] =
false; }
2368 else if(!posSlope && channels[ihit].Cell() < lastHitCell){
2369 trkHits[ihit] =
false;
2372 else if(!posSlope && channels[ihit].Cell() >= lastHitCell){ newTrkHits[ihit] =
false; }
2377 if(nhitsnew > 0){
return; }
2390 for(
int ipl = 0; ipl <
fNChans; ++ipl){
2405 std::vector<bool> cellSortedHits(planeSortedHits.size());
2407 for(
unsigned int i = 0;
i < planeSortedHits.size(); ++
i){
2410 cellSortedHits[cellsort] = planeSortedHits[
i];
2412 planeSortedHits = cellSortedHits;
2418 std::vector<bool> planeSortedHits(cellSortedHits.size());
2420 for(
unsigned int i = 0;
i < cellSortedHits.size(); ++
i){
2429 planeSortedHits[planesort] = cellSortedHits[
i];
2431 cellSortedHits = planeSortedHits;
2437 std::vector<std::array<double, 2>> cellSortedVec(planeSortedVec.size());
2439 for(
unsigned int i = 0;
i < planeSortedVec.size(); ++
i){
2442 cellSortedVec[cellsort] = planeSortedVec[
i];
2444 planeSortedVec = cellSortedVec;
2450 std::vector<std::array<double, 2>> planeSortedVec(cellSortedVec.size());
2452 for(
unsigned int i = 0;
i < cellSortedVec.size(); ++
i){
2461 planeSortedVec[planesort] = cellSortedVec[
i];
2463 cellSortedVec = planeSortedVec;
2471 if(trkHits[
i]){ ++nhits; }
2490 for(
int i = fNChans-1;
i >=0; --
i){
2496 if(minpl == maxpl){
return true; }
2497 else{
return false; }
size_t NTrajectoryPoints() const
int CountMissedCellsOnLine(geo::View_t view, double vi, double zi, double vf, double zf, int maxCells)
Variation of geo::CountCellsOnlineFast function specific to KalmanTrack algorithm.
std::vector< std::vector< double > > avgXPosZRange
void reserve(size_type n)
void PlaneToCellSortHits(std::vector< bool > &planeSortedHits)
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
virtual void SetStart(TVector3 start)
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void GetCenter(double *xyz, double localz=0.0) const
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
void PlaneToCellSort(std::vector< std::array< double, 2 > > &planeSortedHits)
double FilterOnly(const std::vector< bool > &trkHits, std::vector< std::array< double, 2 > > &trkTrajs, std::vector< std::array< double, 2 > > &trkDirs, std::vector< double > &initialP, int zDir, bool zProp, int &newOutlierPos)
A collection of geometry functions used by the KalmanTrack tracking algorithm.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
void SetHit(art::Ptr< rb::CellHit > chit)
std::vector< rb::Track > FindTracks(KalmanGeoHelper &kgeo, std::vector< trk::LocatedChan > sliceChans)
void CheckTrack(KalmanGeoHelper &kgeo, std::vector< bool > &trkHits, std::vector< std::array< double, 2 > > trkTrajs, bool zProp, std::vector< bool > &newTrkHits)
std::string fClusterInput
Vertical planes which measure X.
void reconfigure(const fhicl::ParameterSet &p)
unsigned int Ncells() const
Number of cells in this plane.
KalmanTrack(fhicl::ParameterSet const &pset)
A collection of associated CellHits.
std::vector< double > avgXPos
A rb::Prong with full reconstructed trajectory.
const PlaneGeo * Plane(unsigned int i) const
Representation of a detector channel with positional information.
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
unsigned distance(const T &t1, const T &t2)
Horizontal planes which measure Y.
void SortByCell(std::vector< trk::LocatedChan > &c)
Module that kips a configurable number of events between each that it allows through. Note that this module really skips (N-1) events, it uses a simple modular division as its critera. This module will cut down the data sample to 1/N of its original size.
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
ProductID put(std::unique_ptr< PROD > &&product)
std::vector< trk::LocatedChan > fSliceChansCellSort
#define P(a, b, c, d, e, x)
unsigned short Cell() const
Track finder for cosmic rays.
View_t View() const
Which coordinate does this plane measure.
std::map< ToFCounter, std::vector< unsigned int > > channels
double CalcMissFrac(double dist, double w, geo::View_t view)
Estimation of the mip missed hit fraction as a function of the dist through a cell and cell depth w...
void RemoveHitsFromSignal(const rb::Track &track)
void push_back(Ptr< U > const &p)
double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
void AppendTrajectoryPoint(TVector3 pt)
T get(std::string const &key) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
art::ServiceHandle< chaninfo::BadChanList > fbadc
Geometry information for a single readout plane.
void produce(art::Event &evt)
std::vector< trk::LocatedChan > fSliceChans
double SingleSegment(std::vector< std::vector< double > > &iP, std::vector< bool > &propDir, std::vector< std::vector< bool > > &trkHits, std::vector< std::vector< std::array< double, 2 > > > &trkTraj, std::vector< std::vector< std::array< double, 2 > > > &trkDirs, int firsthit)
void CellToPlaneSort(std::vector< std::array< double, 2 > > &cellSortedHits)
void ClearTrajectoryPoints()
Forget about all trajectory points.
Collect Geo headers and supply basic geometry functions.
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
double DetHalfHeight() const
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
std::vector< double > avgYPos
void RemoveSameTracks(std::vector< std::vector< bool > > &trkHits, std::vector< std::vector< std::array< double, 2 > > > &trkTrajs, std::vector< std::vector< std::array< double, 2 > > > &trkDirs, std::vector< double > &chi2, std::vector< std::vector< double > > &pfiltState, std::vector< bool > &zProp, std::vector< int > &ntrkHits, std::vector< bool > &ignoreTrk)
std::vector< std::vector< double > > avgYPosZRange
int GetLastHit(const std::vector< bool > &trkHits)
std::vector< int > fPlaneToCell
bool IsSinglePlane(const std::vector< bool > &trkHits)
double DetHalfWidth() const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
art::ServiceHandle< geo::Geometry > fgeom
const unsigned int NextPlaneInView(unsigned int p, int d=+1) const
void SortByPlane(std::vector< trk::LocatedChan > &c)
int CountHits(const std::vector< bool > &trkHits)
fvar< T > floor(const fvar< T > &x)
int GetFirstHit(const std::vector< bool > &trkHits)
void RemoveTrajectoryPoint(unsigned int i)
Remove the ith trajectory point from the track.
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
double BestTrackPoint(double z, double v, double halfz, double halfv, std::array< double, 2 > &bestPos, double z0, double v0, double dz, double dv)
Estimate the position on a line within a box.
unsigned int NPlanes() const
TVector3 Stop() const
Position of the final trajectory point.
Encapsulate the cell geometry.
bool IsBad(int plane, int cell)
int SegmentFinder(KalmanGeoHelper &kgeo, std::vector< std::vector< double > > &iP, std::vector< bool > &propDir, std::vector< std::vector< bool > > &trkHits, std::vector< std::vector< std::array< double, 2 > > > &trkTraj, std::vector< std::vector< std::array< double, 2 > > > &trkDirs)
void CellToPlaneSortHits(std::vector< bool > &cellSortedHits)
double FilterTracker(KalmanGeoHelper &kgeo, int zDir, bool zProp, std::vector< bool > &trkHits, std::vector< std::array< double, 2 > > &trkTrajs, std::vector< std::array< double, 2 > > &trkDirs, std::vector< double > &iP, int &outnhits, int firsthit=-1)
bool IsLowEfficiency(int plane, int cell)