1121 if(nHits == 0){
return 0.0; }
1126 if(zDir < 0){ iFirstHit =
fNChans-1; }
1128 bool firstHitfound =
false;
1129 while(!firstHitfound && iFirstHit<fNChans && iFirstHit >= 0){
1131 if(trkHits[iFirstHit]){ firstHitfound =
true; }
1132 else{iFirstHit+=zDir;}
1135 if(trkHits[iFirstHit]){ firstHitfound =
true; }
1136 else{ iFirstHit+=zDir; }
1140 else{ iFirstHit = firsthit; }
1142 unsigned short iPlane = channels[iFirstHit].Plane();
1143 unsigned short iCell = channels[iFirstHit].Cell();;
1146 int iCurrentHit = 0;
1147 if(zDir < 0){ iCurrentHit =
fNChans-1; }
1148 bool currentHitfound =
false;
1149 while(!currentHitfound && iCurrentHit<fNChans && iCurrentHit >= 0){
1151 if(channels[iCurrentHit].Plane() == iPlane){ currentHitfound =
true; }
1152 else{iCurrentHit+=zDir;}
1155 if(channels[iCurrentHit].Cell() == iCell){ currentHitfound =
true; }
1156 else{ iCurrentHit+=zDir; }
1161 const std::array<double, 2> &
start = trkTrajs[iFirstHit];
1162 std::array<double, 2>
dir = trkDirs[iFirstHit];
1163 int iCurrentTrackHit = iFirstHit;
1169 position = start[0];
1170 otherpos = start[1];
1171 if(dir[1] == 0){ dir[1] = 0.0001; }
1172 slope = dir[0]/dir[1];
1175 position = start[1];
1176 otherpos = start[0];
1177 if(dir[0] == 0){ dir[0] = 0.0001; }
1178 slope = dir[1]/dir[0];
1181 double R = cellVarW/3.0;
1182 if(!zProp){ R = cellVarL/3.0; };
1187 int iNextHit = iCurrentHit;
1196 else{ outnhits = nHits; }
1198 int nTrackHits = outnhits;
1208 short currentHitPlane = iPlane;
1209 int expectedHitPlane = currentHitPlane;
1210 short currentHitCell = iCell;
1211 short expectedHitCell = currentHitCell;
1213 std::array<double, 2> bestpos;
1215 while(((iCurrentHit < (
fNChans-1) && zDir > 0) || (iCurrentHit > 0 && zDir < 0)) && !done){
1217 short nextHitPlane = channels[iNextHit].Plane();
1218 short nextHitCell = channels[iNextHit].Cell();
1222 currentHitPlane = channels[iCurrentTrackHit].Plane();
1223 currentHitCell = channels[iCurrentTrackHit].Cell();
1227 if(zProp){sameLayerTest = ((nextHitPlane-currentHitPlane) == (expectedHitPlane-currentHitPlane)); }
1228 else{ sameLayerTest = (((nextHitCell-currentHitCell) == (expectedHitCell-currentHitCell))); }
1233 double peWeight = 0;
1234 double layerState[3] = {state[0], state[1], state[2]};
1235 double layerP[3] = {errorCov[0], errorCov[1], errorCov[2]};
1238 int firstLayerHit = iNextHit;
1240 if(zProp){
while(k < fNChans && k >= 0 && channels[k].Plane() == nextHitPlane){ k+=zDir; } }
1241 else{
while(k < fNChans && k >= 0 && channels[k].Cell() == nextHitCell){ k+=zDir; } }
1242 int nadd = zDir*(k-iNextHit);
1243 bool layerTrackHits[nadd];
1244 int iLayerHit = firstLayerHit;
1245 int fLayerHit = k-zDir;
1246 for(
int i = 0;
i < nadd; ++
i){
1247 layerTrackHits[
i] =
false;
1248 if(trkHits[iNextHit]){
1249 trkHits[iNextHit] =
false;
1252 iCurrentHit = iNextHit;
1253 iNextHit = iCurrentHit + zDir;
1255 bool doneAdding =
false;
1256 if(nadd == 0){ doneAdding =
true;}
1259 int layerWindowlow = 0;
1260 int layerWindowhigh = 999999;
1262 if(((state[2] <= 0 && zDir < 0) || (state[2] >= 0 && zDir > 0))){
1264 layerWindowlow = channels[iCurrentTrackHit].Plane();
1265 if(
abs(nextHitCell-currentHitCell)>1){
1266 layerWindowlow = channels[iCurrentTrackHit].Plane()-2;
1267 layerWindowhigh = channels[iCurrentTrackHit].Plane()+6;
1273 layerWindowhigh = channels[iCurrentTrackHit].Plane();
1274 if(
abs(nextHitCell-currentHitCell)>1){
1275 layerWindowhigh = channels[iCurrentTrackHit].Plane()+2;
1276 layerWindowlow = channels[iCurrentTrackHit].Plane()-6;
1282 double minDistBest = -1;
1283 int bestHitIndex = 9999;
1284 if(nadd == 1 && !layerTrackHits[0]){ bestHitIndex = 0; }
1286 for(
int i = 0;
i<nadd;++
i){
1288 if(!layerTrackHits[
i]){
1289 int currChanIdx = firstLayerHit+i*zDir;
1292 if(zProp && (channels[currChanIdx].Cell()<layerWindowlow || channels[currChanIdx].Cell()>layerWindowhigh)){
continue; }
1293 else if(!zProp && (channels[currChanIdx].Plane()<layerWindowlow || channels[currChanIdx].Plane()>layerWindowhigh)){
continue; }
1297 channels[currChanIdx].GetCenter(xyzHit);
1298 if(zProp){kgeo.BestTrackPoint(xyzHit[2],xyzHit[
fView],fCellL,fCellW,bestpos,layerState[1],layerState[0],dir[1],dir[0]);}
1299 else{kgeo.BestTrackPoint(xyzHit[2],xyzHit[
fView],fCellL,fCellW,bestpos,layerState[0],layerState[1],dir[1],dir[0]);}
1301 double minHitDist = 0.0;
1302 if(zProp){ minHitDist =
abs(bestpos[1]-(layerState[0]+layerState[2]*(bestpos[0]-layerState[1])));}
1303 else{ minHitDist =
abs(bestpos[0]-(layerState[0]+layerState[2]*(bestpos[1]-layerState[1]))); }
1304 if((minHitDist < minDistBest) || minDistBest < 0){
1305 bestHitIndex = i*zDir;
1306 minDistBest = minHitDist;
1311 if(bestHitIndex != 9999){
1312 int bestChanIdx = firstLayerHit+bestHitIndex;
1314 channels[bestChanIdx].GetCenter(xyzHit);
1317 if(zProp){ delta = xyzHit[2] - layerState[1]; }
1318 else{ delta = xyzHit[
fView] - layerState[1]; }
1319 Q[0] = delta*delta*Q[2] + layerP[0];
1323 xMinus[0] = layerState[0] + delta*layerState[2];
1324 xMinus[2] = layerState[2];
1325 pMinus[0] = layerP[0] + 2.0*delta*layerP[1] + delta*delta*layerP[2] + Q[0];
1326 pMinus[1] = layerP[1] + delta*layerP[2] + Q[1];
1327 pMinus[2] = layerP[2] + Q[2];
1330 double denom = (pMinus[0]+
R);
1331 K[0] = pMinus[0]/denom;
1332 K[1] = pMinus[1]/denom;
1335 double candidatePos;
1338 candidatePos = xyzHit[
fView];
1343 candidatePos = xyzHit[2];
1344 ipos = xyzHit[
fView];
1349 xPredict[0] = xMinus[0]+K[0]*(candidatePos-xMinus[0]);
1350 xPredict[2] = xMinus[2]+K[1]*(candidatePos-xMinus[0]);
1351 cPredict[0] = pMinus[0]*R/denom;
1352 cPredict[1] = pMinus[1]*R/denom;
1353 cPredict[2] = (pMinus[2]*pMinus[0]-pMinus[1]*pMinus[1]+pMinus[2]*
R)/denom;
1355 double r2 = (candidatePos - xPredict[0])*(candidatePos - xPredict[0]);
1356 double sl2 = xPredict[2]*xPredict[2];
1357 double deltaChi = 0.0;
1359 deltaChi = (r2*(1+sl2)*(1+sl2))/
1360 (((1+sl2)*(1+sl2)*(cellVarW/3+sl2*cellVarL/3+cPredict[0]))+r2*sl2*cPredict[2]/2);
1363 deltaChi = (r2*(1+sl2)*(1+sl2))/
1364 (((1+sl2)*(1+sl2)*(cellVarL/3+sl2*cellVarW/3+cPredict[0]))+r2*sl2*cPredict[2]/2);
1367 if(deltaChi <= maxDeltaChi){
1370 trkHits[bestChanIdx] =
true;
1373 layerTrackHits[(bestHitIndex/zDir)] =
true;
1374 currentHitPlane = channels[bestChanIdx].Plane();
1375 currentHitCell = channels[bestChanIdx].Cell();
1376 float currentPE = channels[bestChanIdx].GetHit()->PE();
1377 iCurrentTrackHit = bestChanIdx;
1379 layerState[0] = xPredict[0];
1380 layerState[1] = xPredict[1];
1381 layerState[2] = xPredict[2];
1382 layerP[0] = cPredict[0];
1383 layerP[1] = cPredict[1];
1384 layerP[2] = cPredict[2];
1385 peWeight = currentPE;
1388 layerState[0] = (layerState[0]*peWeight+xPredict[0]*currentPE)/(peWeight+currentPE);
1389 layerState[1] = (layerState[1]*peWeight+xPredict[1]*currentPE)/(peWeight+currentPE);
1390 layerState[2] = (layerState[2]*peWeight+xPredict[2]*currentPE)/(peWeight+currentPE);
1391 layerP[0] = (layerP[0]*peWeight+cPredict[0]*currentPE)/(peWeight+currentPE);
1392 layerP[1] = (layerP[1]*peWeight+cPredict[1]*currentPE)/(peWeight+currentPE);
1393 layerP[2] = (layerP[2]*peWeight+cPredict[2]*currentPE)/(peWeight+currentPE);
1394 peWeight+=currentPE;
1398 dir[0] = layerState[2]/
sqrt(1.0+layerState[2]*layerState[2]);
1399 dir[1] = 1.0/
sqrt(1.0+layerState[2]*layerState[2]);
1401 int newLayerWindowlow = currentHitCell-
fGapAllowed;
1402 int newLayerWindowhigh = currentHitCell+
fGapAllowed;
1404 while(newLayerWindowlow > 0 &&
1405 (
fbadc->
IsBad(currentHitPlane,newLayerWindowlow)
1407 --newLayerWindowlow;
1410 (
fbadc->
IsBad(currentHitPlane,newLayerWindowhigh)
1412 ++newLayerWindowhigh;
1415 layerWindowlow = newLayerWindowlow;
1416 layerWindowhigh = newLayerWindowhigh;
1419 if((currentHitCell - layerWindowlow) < (layerWindowhigh - currentHitCell)){
1423 int newLayerWindowlow = currentHitCell-
fGapAllowed;
1424 while(newLayerWindowlow > 0 &&
1425 (
fbadc->
IsBad(currentHitPlane,newLayerWindowlow)
1427 --newLayerWindowlow;
1429 layerWindowlow = newLayerWindowlow;
1436 int newLayerWindowhigh = currentHitCell+
fGapAllowed;
1438 (
fbadc->
IsBad(currentHitPlane,newLayerWindowhigh)
1440 ++newLayerWindowhigh;
1442 layerWindowhigh = newLayerWindowhigh;
1449 dir[0] = 1.0/
sqrt(1.0+layerState[2]*layerState[2]);
1450 dir[1] = layerState[2]/
sqrt(1.0+layerState[2]*layerState[2]);
1452 if(((state[2] <= 0 && zDir < 0) || (state[2] >= 0 && zDir > 0))){
1453 int newLayerWindowlow = currentHitPlane;
1455 layerWindowlow = newLayerWindowlow;
1459 int newlayerWindowhightest = newLayerWindowhigh;
1461 (
fbadc->
IsBad(newlayerWindowhightest,currentHitCell)
1468 layerWindowhigh = newLayerWindowhigh;
1472 int newLayerWindowhigh = currentHitPlane;
1473 layerWindowhigh = newLayerWindowhigh;
1477 int newlayerWindowlowtest = newLayerWindowlow;
1479 (
fbadc->
IsBad(newlayerWindowlowtest,currentHitCell)
1486 layerWindowlow = newLayerWindowlow;
1490 if(((state[2] <= 0 && zDir < 0) || (state[2] >= 0 && zDir > 0))){
1491 layerWindowlow = currentHitPlane;
1500 int newLayerWindowhightest = newLayerWindowhigh;
1502 (
fbadc->
IsBad(newLayerWindowhightest,currentHitCell)
1508 layerWindowhigh = newLayerWindowhigh;
1512 layerWindowhigh = currentHitPlane;
1521 int newLayerWindowlowtest = newLayerWindowlow;
1523 (
fbadc->
IsBad(newLayerWindowlowtest,currentHitCell)
1529 layerWindowlow = newLayerWindowlow;
1535 else{ layerTrackHits[(bestHitIndex/zDir)] =
true; }
1538 else{ doneAdding =
true;}
1544 state[0] = layerState[0];
1545 state[1] = layerState[1];
1546 state[2] = layerState[2];
1547 errorCov[0] = layerP[0];
1548 errorCov[1] = layerP[1];
1549 errorCov[2] = layerP[2];
1553 if(fLayerHit < iLayerHit){
std::swap(iLayerHit,fLayerHit); }
1554 for(
int ihit = iLayerHit; ihit<=fLayerHit; ++ihit){
1556 trkTrajs[ihit][0] = state[0];
1557 trkTrajs[ihit][1] = state[1];
1558 trkDirs[ihit][0] = dir[0];
1559 trkDirs[ihit][1] = dir[1];
1566 if(fLayerHit < iLayerHit){
std::swap(iLayerHit,fLayerHit); }
1567 for(
int ihit = iLayerHit; ihit<=fLayerHit; ++ihit){
1569 trkTrajs[ihit][0] = state[1];
1570 trkTrajs[ihit][1] = state[0];
1571 trkDirs[ihit][0] = dir[0];
1572 trkDirs[ihit][1] = dir[1];
1575 expectedHitCell = currentHitCell+zDir;
1576 if(expectedHitCell < 0 || expectedHitCell >
NCells[
fView] ){ done =
true; }
1578 if(iNextHit > (
fNChans-1) || iNextHit < 0){ done =
true; }
1584 if(iNextHit > (
fNChans-1) || iNextHit < 0){ done =
true; }
1586 short nextHitPlane = channels[iNextHit].Plane();
1587 short nextHitCell = channels[iNextHit].Cell();
1592 channels[iNextHit].GetCenter(nextpt);
1595 double offset = zDir*(fCellL+2.0);
1596 currtraj[1] = state[1] +
offset;
1597 currtraj[0] = state[2]*offset + state[0];
1600 nextpt[
fView] = state[2]*(nextpt[2]-currtraj[1])+currtraj[0];
1607 double offset = zDir*(fCellW+2.0);
1608 currtraj[0] = state[1] +
offset;
1609 currtraj[1] = state[2]*offset + state[0];
1612 nextpt[2] = state[2]*(nextpt[
fView]-currtraj[0])+currtraj[1];
1613 if(nextpt[2] < 0 || nextpt[2] >
fgeom->
DetLength()){ done =
true; }
1615 misses+=kgeo.CountMissedCellsOnLine(fView,currtraj[0],currtraj[1],nextpt[fView],nextpt[2],fGapAllowed);
1616 if(misses > maxMisses){done =
true;}
1618 expectedHitCell = nextHitCell;
1619 expectedHitPlane = nextHitPlane;
1629 channels[iNextHit].GetCenter(nextpt);
1632 double offset = zDir*(fCellL+2.0);
1633 currtraj[1] = state[1] +
offset;
1634 currtraj[0] = state[2]*offset + state[0];
1637 nextpt[
fView] = state[2]*(nextpt[2]-currtraj[1])+currtraj[0];
1644 double offset = zDir*(fCellW+2.0);
1645 currtraj[0] = state[1] +
offset;
1646 currtraj[1] = state[2]*offset + state[0];
1649 nextpt[2] = state[2]*(nextpt[
fView]-currtraj[0])+currtraj[1];
1650 if(nextpt[2] < 0 || nextpt[2] >
fgeom->
DetLength()){ done =
true; }
1652 misses+=kgeo.CountMissedCellsOnLine(fView,currtraj[0],currtraj[1],nextpt[fView],nextpt[2],fGapAllowed);
1653 if(misses > maxMisses){done =
true;}
1655 expectedHitPlane = nextHitPlane;
1656 expectedHitCell = nextHitCell;
std::vector< trk::LocatedChan > fSliceChansCellSort
std::map< ToFCounter, std::vector< unsigned int > > channels
art::ServiceHandle< chaninfo::BadChanList > fbadc
std::vector< trk::LocatedChan > fSliceChans
double DetHalfHeight() const
double DetHalfWidth() const
art::ServiceHandle< geo::Geometry > fgeom
const unsigned int NextPlaneInView(unsigned int p, int d=+1) const
int CountHits(const std::vector< bool > &trkHits)
bool IsBad(int plane, int cell)
bool IsLowEfficiency(int plane, int cell)