Public Types | Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes | List of all members
trk::KalmanTrack Class Reference
Inheritance diagram for trk::KalmanTrack:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 KalmanTrack (fhicl::ParameterSet const &pset)
 
virtual ~KalmanTrack ()
 
void produce (art::Event &evt)
 
void reconfigure (const fhicl::ParameterSet &p)
 
std::string workerType () const
 
void doBeginJob ()
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Protected Member Functions

ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Private Member Functions

std::vector< rb::Track > FindTracks (KalmanGeoHelper &kgeo, std::vector< trk::LocatedChan > sliceChans)
 
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)
 
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)
 
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)
 
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)
 
void RemoveHitsFromSignal (const rb::Track &track)
 
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)
 
void CheckTrack (KalmanGeoHelper &kgeo, std::vector< bool > &trkHits, std::vector< std::array< double, 2 > > trkTrajs, bool zProp, std::vector< bool > &newTrkHits)
 
void PlaneToCellSortHits (std::vector< bool > &planeSortedHits)
 
void CellToPlaneSortHits (std::vector< bool > &cellSortedHits)
 
void PlaneToCellSort (std::vector< std::array< double, 2 > > &planeSortedHits)
 
void CellToPlaneSort (std::vector< std::array< double, 2 > > &cellSortedHits)
 
void CreateHitMaps ()
 
int CountHits (const std::vector< bool > &trkHits)
 
bool IsSinglePlane (const std::vector< bool > &trkHits)
 
int GetFirstHit (const std::vector< bool > &trkHits)
 
int GetLastHit (const std::vector< bool > &trkHits)
 

Private Attributes

double fMaxSliceHits
 
double fMaxHitCut
 
int fWhoaNelly
 
std::string fClusterInput
 
double fDeltaChiAllowed
 
unsigned int fMinHits
 
int fGapAllowed
 
int fLargeSliceHits
 
bool fBadChannels
 
bool fLongTrackLoop
 
double fMinGapProb
 
bool fObeyPreselection
 
art::ServiceHandle< geo::Geometry > fgeom
 
art::ServiceHandle< chaninfo::BadChanList > fbadc
 
std::vector< double > avgXPos
 
std::vector< std::vector< double > > avgXPosZRange
 
std::vector< double > avgYPos
 
std::vector< std::vector< double > > avgYPosZRange
 
double avgX
 
double avgY
 
double fSysVariance = 0.0001
 
int NCells [2]
 
int fNChans
 
std::vector< trk::LocatedChan > fSliceChans
 
std::vector< trk::LocatedChan > fSliceChansCellSort
 
std::vector< int > fPlaneToCell
 
geo::View_t fView
 
double fCellW
 
double fCellL
 
bool fTrySingle
 

Detailed Description

Definition at line 32 of file KalmanTrack_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 17 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 25 of file Producer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 18 of file EDProducer.h.

Constructor & Destructor Documentation

trk::KalmanTrack::KalmanTrack ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 147 of file KalmanTrack_module.cc.

References reconfigure().

148  : EDProducer(pset)
149  {
150  reconfigure(pset);
151  produces< std::vector<rb::Track> >();
152  produces< art::Assns<rb::Track, rb::Cluster> >();
153  }
art::EDProducer::EDProducer
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
trk::KalmanTrack::reconfigure
void reconfigure(const fhicl::ParameterSet &p)
Definition: KalmanTrack_module.cc:162
trk::KalmanTrack::~KalmanTrack ( )
virtual

Definition at line 156 of file KalmanTrack_module.cc.

157  {
158  }

Member Function Documentation

void trk::KalmanTrack::CellToPlaneSort ( std::vector< std::array< double, 2 > > &  cellSortedHits)
private

Definition at line 2449 of file KalmanTrack_module.cc.

References fPlaneToCell, MECModelEnuComparisons::i, and calib::j.

Referenced by FindTracks().

2450  {
2451  std::vector<std::array<double, 2>> planeSortedVec(cellSortedVec.size());
2452  // loop over the cellSortedVec and set the corresponding planeSortedHit to true
2453  for(unsigned int i = 0; i < cellSortedVec.size(); ++i){
2454  int planesort = 0;
2455  // get the corresponding plane sorted hit location
2456  for(unsigned int j = 0; j < fPlaneToCell.size(); ++j){
2457  if(fPlaneToCell[j] == (int)i){
2458  planesort = j;
2459  break;
2460  }
2461  }
2462  planeSortedVec[planesort] = cellSortedVec[i];
2463  }
2464  cellSortedVec = planeSortedVec;
2465  }
calib::j
const double j
Definition: BetheBloch.cxx:29
trk::KalmanTrack::fPlaneToCell
std::vector< int > fPlaneToCell
Definition: KalmanTrack_module.cc:67
void trk::KalmanTrack::CellToPlaneSortHits ( std::vector< bool > &  cellSortedHits)
private

Definition at line 2417 of file KalmanTrack_module.cc.

References fPlaneToCell, MECModelEnuComparisons::i, and calib::j.

Referenced by FindTracks().

2418  {
2419  std::vector<bool> planeSortedHits(cellSortedHits.size());
2420  // loop over the cellSortedHits and set the corresponding planeSortedHit to true
2421  for(unsigned int i = 0; i < cellSortedHits.size(); ++i){
2422  int planesort = 0;
2423  // get the corresponding plane sorted hit location
2424  for(unsigned int j = 0; j < fPlaneToCell.size(); ++j){
2425  if(fPlaneToCell[j] == (int)i){
2426  planesort = j;
2427  break;
2428  }
2429  }
2430  planeSortedHits[planesort] = cellSortedHits[i];
2431  }
2432  cellSortedHits = planeSortedHits;
2433  }
calib::j
const double j
Definition: BetheBloch.cxx:29
trk::KalmanTrack::fPlaneToCell
std::vector< int > fPlaneToCell
Definition: KalmanTrack_module.cc:67
void trk::KalmanTrack::CheckTrack ( KalmanGeoHelper &  kgeo,
std::vector< bool > &  trkHits,
std::vector< std::array< double, 2 > >  trkTrajs,
bool  zProp,
std::vector< bool > &  newTrkHits 
)
private

Definition at line 2116 of file KalmanTrack_module.cc.

References avgXPos, avgXPosZRange, avgYPos, avgYPosZRange, trk::KalmanGeoHelper::CalcMissFrac(), channels, CountHits(), trk::KalmanGeoHelper::CountMissedCellsOnLine(), distance(), fbadc, fBadChannels, fMinGapProb, fNChans, fSliceChans, fSliceChansCellSort, fView, MECModelEnuComparisons::i, makeTrainCVSamples::int, chaninfo::BadChanList::IsBad(), chaninfo::BadChanList::IsLowEfficiency(), geo::kX, geo::kY, NDAPDHVSetting::plane, and w.

Referenced by FindTracks().

2119  {
2120  newTrkHits.resize(fNChans);
2121  std::vector<trk::LocatedChan> & channels = zProp? fSliceChans : fSliceChansCellSort;
2122  // step through this track and make sure that there aren't any gaps larger than allowed
2123  // If gaps exist break track into two seperate pieces
2124  double minProb = fMinGapProb;
2125 
2126  // first get a list of hits that the track has in it
2127  int nhits = CountHits(trkHits);
2128 
2129  std::vector<std::array<double, 2> > trjPointsNew(nhits);
2130  std::vector<int> trkHitsIdx(nhits);
2131  int ihit = 0;
2132  for(int i = 0; i < fNChans; ++i){
2133  if(trkHits[i]){
2134  trjPointsNew[ihit] = trkTrajs[i];
2135  trkHitsIdx[ihit] = i;
2136  ++ihit;
2137  }
2138  }
2139 
2140  unsigned int lastPlane = 10000; // last plane of a missed hit
2141  unsigned int lastCell = 10000; // last cell of a missed hit
2142  unsigned int lastHitPlane = 10000; // last plane with a hit in it
2143  unsigned int lastHitCell = 10000; // last cell with a hit in it
2144  bool posSlope = true; // does this track have a positive slope, need this for sorting
2145  bool spillover = false; // does a gap exist between last trajectory point to this point that needs to be added
2146  int spillovergap = 0; // last gap in cells that spills over into current gap
2147  double spilloverprob = 1.0; // last gap in probability that spills over into current gap
2148  int allMissed = 0; // total number of missed cells
2149  int biggestGap = 0; // biggest gap in the track
2150 
2151  // loop over all trajectory points
2152  for(int itraj = 0; itraj < (nhits-1); ++itraj){
2153 
2154  // check if this trajectory point is the same as the next, can't possibly have gap in these cases
2155  if(trjPointsNew[itraj][0] == trjPointsNew[itraj+1][0] && trjPointsNew[itraj][1] == trjPointsNew[itraj+1][1]){ continue; }
2156 
2157  bool sort = false;
2158  std::vector<geo::OfflineChan> missedCells; // missed cells
2159  std::vector<double> distance; // 2d distance in missed cells
2160 
2161  // total missed cells between this and next trajectory point
2162  int totMissed = kgeo.CountMissedCellsOnLine(fView,trjPointsNew[itraj][0],trjPointsNew[itraj][1],
2163  trjPointsNew[itraj+1][0],trjPointsNew[itraj+1][1],missedCells,distance);
2164  // If didn't go through any cells, no gaps to check and nothing to update
2165  if(totMissed == 0){ continue; }
2166 
2167  // order the hits in the order they appear along the track path, only neccessary if slope is negative
2168  if(trjPointsNew[itraj][1] != trjPointsNew[itraj+1][1] &&
2169  ((trjPointsNew[itraj+1][0]-trjPointsNew[itraj][0])/(trjPointsNew[itraj+1][1]-trjPointsNew[itraj][1])) < 0){ sort = true; }
2170 
2171  if(sort){
2172  // already ordered by plane, now do it by cell within the plane
2173  posSlope = false;
2174  std::vector<geo::OfflineChan> orderedMissedCells(totMissed);
2175  std::vector<double> orderedDistance(totMissed);
2176  unsigned int firstCellIndex = 0;
2177  unsigned int lastCellIndex = 0;
2178  unsigned int fillCellIndex = 0;
2179  unsigned int plane = missedCells[0].Plane();
2180  for(int iCell = 0; iCell<totMissed;++iCell){
2181  if(missedCells[iCell].Plane() != (int)plane || (int)iCell == totMissed-1){
2182  if(iCell == totMissed-1 && missedCells[iCell].Plane() == (int)plane){lastCellIndex = iCell;}
2183  // first fill the ordered cells vector
2184  for(int fillCell = (int) lastCellIndex; fillCell>= (int) firstCellIndex;--fillCell){
2185  orderedMissedCells[fillCellIndex] = missedCells[fillCell];
2186  orderedDistance[fillCellIndex] = distance[fillCell];
2187  ++fillCellIndex;
2188  }
2189  if(iCell == totMissed-1 && missedCells[iCell].Plane() != (int)plane){
2190  orderedMissedCells[totMissed-1] = missedCells[totMissed-1];
2191  orderedDistance[totMissed-1] = distance[totMissed-1];
2192  lastCellIndex = iCell;
2193  if(missedCells[iCell].Plane() != plane){firstCellIndex = iCell;}
2194  }
2195  // update the counters and plane
2196  firstCellIndex = iCell;
2197  lastCellIndex = iCell;
2198  plane = missedCells[iCell].Plane();
2199  }
2200  else{ lastCellIndex = iCell; }
2201  }
2202  missedCells = orderedMissedCells;
2203  distance = orderedDistance;
2204  }
2205 
2206  int currentgap = 0; // current gap in cells between this trajectory point and the next
2207  double currprob = 1.0; // current probability of cell gap between this and next trajectory point
2208 
2209  // if there is a spillover add it to the current gap
2210  if(spillover){
2211  currentgap+=spillovergap;
2212  currprob = spilloverprob;
2213  }
2214 
2215  // now the cells are ordered so we can look for gaps
2216  for(int iCell = 0; iCell < totMissed; ++iCell){
2217 
2218  bool ontrack = false; // Is this cell a hit on the track?
2219  bool inbadc = false; // Is this cell a bad channel?
2220 
2221  // check if the track has a hit in this cell
2222  for(int ihit = 0; ihit < nhits; ++ihit){
2223  if(channels[trkHitsIdx[ihit]].Plane() == missedCells[iCell].Plane() &&
2224  channels[trkHitsIdx[ihit]].Cell() == missedCells[iCell].Cell()){
2225  ontrack = true;
2226  lastHitPlane = missedCells[iCell].Plane();
2227  lastHitCell = missedCells[iCell].Cell();
2228  break;
2229  }
2230  }
2231 
2232  // check if this cell was a bad channel, if so don't count against the track, but also don't reset gaps
2233  if(fBadChannels && !ontrack){
2234  if(fbadc->IsBad(missedCells[iCell].Plane(),missedCells[iCell].Cell()) ||
2235  fbadc->IsLowEfficiency(missedCells[iCell].Plane(),missedCells[iCell].Cell()) ){ inbadc = true; }
2236  }
2237 
2238  // if hit is on track and not in a bad channel reset current and spillover gaps to zero
2239  if(ontrack){
2240  currentgap = 0;
2241  currprob = 1.0;
2242  spillovergap = 0;
2243  spilloverprob = 1.0;
2244  spillover = false;
2245  }
2246  else if(!ontrack && !inbadc){
2247  // if cell is not on the track and its not a bad channel calculate the gap probabilities
2248  double lenInCell = distance[iCell];
2249  double zave = 0.5*trjPointsNew[itraj][1]+0.5*trjPointsNew[itraj+1][1];
2250 
2251  // find the approximate w of this cell
2252  double w = 0;
2253  if(fView == geo::kX){
2254  // loop over all of the breakpoints and find what w to use
2255  bool foundrng = false;
2256  if(zave < avgYPosZRange[0][0]){
2257  // use the first range
2258  w = avgYPos[0];
2259  foundrng = true;
2260  }
2261  else if(zave > avgYPosZRange[avgYPosZRange.size()-1][1]){
2262  // use the last range
2263  w = avgYPos[avgYPos.size()-1];
2264  foundrng = true;
2265  }
2266  else{
2267  // loop over all the ranges and find the one that encloses this z
2268  for(unsigned int i = 0; i < avgYPosZRange.size();++i){
2269  if(zave >= avgYPosZRange[i][0] && i <= avgYPosZRange[i][1]){
2270  w = avgYPos[i];
2271  foundrng = true;
2272  break;
2273  }
2274  }
2275  }
2276  if(!foundrng){ w = avgYPos[0]; }
2277  }
2278  else if(fView == geo::kY){
2279  // loop over all of the breakpoints and find what w to use
2280  bool foundrng = false;
2281  if(zave < avgXPosZRange[0][0]){
2282  // use the first range
2283  w = avgXPos[0];
2284  foundrng = true;
2285  }
2286  else if(zave > avgXPosZRange[avgXPosZRange.size()-1][1]){
2287  // use the last range
2288  w = avgXPos[avgXPos.size()-1];
2289  foundrng = true;
2290  }
2291  else{
2292  // loop over all the ranges and find the one that encloses this z
2293  for(unsigned int i = 0; i < avgXPosZRange.size();++i){
2294  if(zave >= avgXPosZRange[i][0] && i<= avgXPosZRange[i][1]){
2295  w = avgXPos[i];
2296  foundrng = true;
2297  break;
2298  }
2299  }
2300  }
2301  if(!foundrng){ w = avgXPos[0]; }
2302  }
2303 
2304  double probMiss = kgeo.CalcMissFrac(lenInCell,w,fView);
2305 
2306  if(probMiss < 1.0){
2307  if(!spillover){
2308  ++currentgap;
2309  currprob*=probMiss;
2310  ++allMissed;
2311  }
2312  else{
2313  // Only add if this doesn't correspond to the last plane/cell of the spillover (don't double count)
2314  if(missedCells[iCell].Plane() != lastPlane ||
2315  missedCells[iCell].Cell() != lastCell){
2316  ++currentgap;
2317  currprob*=probMiss;
2318  ++allMissed;
2319  }
2320  }
2321  }
2322 
2323  //update the gap that spills over into the next track section if it is the last cell
2324  if(iCell == (totMissed-1)){
2325  spillovergap = currentgap;
2326  spilloverprob = currprob;
2327  lastPlane = missedCells[iCell].Plane();
2328  lastCell = missedCells[iCell].Cell();
2329  spillover = true;
2330  }
2331 
2332  } // end of if(!ontrack && !inbadc)
2333 
2334  // if we have crossed the threshold of the gap probability stop, break the track here
2335  if(currprob < minProb && currentgap > 1){
2336  if(lastHitPlane == 10000){
2337  lastHitPlane = missedCells[iCell].Plane();
2338  lastHitCell = missedCells[iCell].Cell();
2339  }
2340  break;
2341  }
2342 
2343  } // end of loop over iCell
2344 
2345  if(currentgap > biggestGap){ biggestGap = currentgap; }
2346 
2347  if(currprob < minProb && currentgap > 1){
2348 
2349  // find all the hits to take off this track
2350  newTrkHits = trkHits;
2351  int nhitsnew = 0;
2352  for(int ihit = 0; ihit < fNChans; ++ihit){
2353  // don't need to do anything if this hit isn't on the track
2354  if(!trkHits[ihit]){ continue; }
2355  // if hit after the last hit plane, it shouldn't be anymore. Add it to new track hits
2356  if(channels[ihit].Plane() > lastHitPlane){
2357  trkHits[ihit] = false;
2358  ++nhitsnew;
2359  }
2360  // if hit is before last hit plane, it should only be on the old track
2361  if(channels[ihit].Plane() < lastHitPlane){ newTrkHits[ihit] = false; }
2362  // if hit is on the last hit plane it should only be on the new track if after the last hit cell
2363  if(channels[ihit].Plane() == lastHitPlane){
2364  if(posSlope && channels[ihit].Cell() > lastHitCell){
2365  trkHits[ihit] = false;
2366  ++nhitsnew;
2367  }
2368  else if(posSlope && channels[ihit].Cell() <= lastHitCell){ newTrkHits[ihit] = false; }
2369  else if(!posSlope && channels[ihit].Cell() < lastHitCell){
2370  trkHits[ihit] = false;
2371  ++nhitsnew;
2372  }
2373  else if(!posSlope && channels[ihit].Cell() >= lastHitCell){ newTrkHits[ihit] = false; }
2374  }
2375  }
2376 
2377  // stop after a gap is found
2378  if(nhitsnew > 0){ return; }
2379  }
2380 
2381  } // end of loop over trajectory points
2382 
2383 
2384  } // end of CheckTracks function
trk::KalmanTrack::avgXPosZRange
std::vector< std::vector< double > > avgXPosZRange
Definition: KalmanTrack_module.cc:57
geo::kX
Vertical planes which measure X.
Definition: PlaneGeo.h:28
trk::KalmanTrack::avgXPos
std::vector< double > avgXPos
Definition: KalmanTrack_module.cc:56
distance
unsigned distance(const T &t1, const T &t2)
Definition: Multiplet_module.cc:33
geo::kY
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
trk::KalmanTrack::fSliceChansCellSort
std::vector< trk::LocatedChan > fSliceChansCellSort
Definition: KalmanTrack_module.cc:66
channels
std::map< ToFCounter, std::vector< unsigned int > > channels
Definition: ToFPositionRecoAnalysis_module.cc:81
trk::KalmanTrack::fbadc
art::ServiceHandle< chaninfo::BadChanList > fbadc
Definition: KalmanTrack_module.cc:55
trk::KalmanTrack::fSliceChans
std::vector< trk::LocatedChan > fSliceChans
Definition: KalmanTrack_module.cc:65
trk::KalmanTrack::avgYPos
std::vector< double > avgYPos
Definition: KalmanTrack_module.cc:58
trk::KalmanTrack::avgYPosZRange
std::vector< std::vector< double > > avgYPosZRange
Definition: KalmanTrack_module.cc:59
trk::KalmanTrack::CountHits
int CountHits(const std::vector< bool > &trkHits)
Definition: KalmanTrack_module.cc:2468
chaninfo::BadChanList::IsBad
bool IsBad(int plane, int cell)
Definition: BadChanList_service.cc:563
w
Float_t w
Definition: plot.C:20
chaninfo::BadChanList::IsLowEfficiency
bool IsLowEfficiency(int plane, int cell)
Definition: BadChanList_service.cc:691
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::consumes ( InputTag const &  tag)
protectedinherited

Definition at line 55 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::consumes(), T, and getGoodRuns4SAM::tag.

56  {
57  return collector_.consumes<T, BT>(tag);
58  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:50
art::ConsumesCollector::consumes
ProductToken< T > consumes(InputTag const &)
Definition: ConsumesCollector.h:44
T
double T
Definition: Xdiff_gwt.C:5
ConsumesCollector& art::ModuleBase::consumesCollector ( )
protectedinherited
template<typename T , BranchType BT>
void art::ModuleBase::consumesMany ( )
protectedinherited

Definition at line 69 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::consumesMany(), and T.

70  {
71  collector_.consumesMany<T, BT>();
72  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:50
T
double T
Definition: Xdiff_gwt.C:5
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::consumesView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::consumesView ( InputTag const &  tag)
inherited

Definition at line 62 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::consumesView(), T, and getGoodRuns4SAM::tag.

63  {
64  return collector_.consumesView<T, BT>(tag);
65  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:50
art::ConsumesCollector::consumesView
ViewToken< Element > consumesView(InputTag const &)
T
double T
Definition: Xdiff_gwt.C:5
int trk::KalmanTrack::CountHits ( const std::vector< bool > &  trkHits)
private

Definition at line 2468 of file KalmanTrack_module.cc.

References fNChans, and MECModelEnuComparisons::i.

Referenced by CheckTrack(), FilterOnly(), FilterTracker(), and FindTracks().

2469  {
2470  int nhits = 0;
2471  for(int i = 0; i < fNChans; ++i){
2472  if(trkHits[i]){ ++nhits; }
2473  }
2474  return nhits;
2475  }
void trk::KalmanTrack::CreateHitMaps ( )
private

Definition at line 2387 of file KalmanTrack_module.cc.

References fNChans, fPlaneToCell, fSliceChans, and fSliceChansCellSort.

Referenced by FindTracks().

2388  {
2389  fPlaneToCell.resize(fNChans);
2390  // loop over all the plane sorted hits
2391  for(int ipl = 0; ipl < fNChans; ++ipl){
2392  art::Ptr<rb::CellHit> plhit = fSliceChans[ipl].GetHit();
2393  // loop over all the cell sorted hits and compare to this hit
2394  for(unsigned int icl = 0; icl < fSliceChansCellSort.size(); ++icl){
2395  if(plhit == fSliceChansCellSort[icl].GetHit()){
2396  fPlaneToCell[ipl] = icl;
2397  break;
2398  }
2399  }
2400  }
2401  } // end of CreateHitMaps()
trk::KalmanTrack::fSliceChansCellSort
std::vector< trk::LocatedChan > fSliceChansCellSort
Definition: KalmanTrack_module.cc:66
trk::KalmanTrack::fSliceChans
std::vector< trk::LocatedChan > fSliceChans
Definition: KalmanTrack_module.cc:65
trk::KalmanTrack::fPlaneToCell
std::vector< int > fPlaneToCell
Definition: KalmanTrack_module.cc:67
void art::detail::Producer::doBeginJob ( )
inherited
bool art::detail::Producer::doBeginRun ( RunPrincipal &  rp,
ModuleContext const &  mc 
)
inherited
bool art::detail::Producer::doBeginSubRun ( SubRunPrincipal &  srp,
ModuleContext const &  mc 
)
inherited
void art::detail::Producer::doEndJob ( )
inherited
bool art::detail::Producer::doEndRun ( RunPrincipal &  rp,
ModuleContext const &  mc 
)
inherited
bool art::detail::Producer::doEndSubRun ( SubRunPrincipal &  srp,
ModuleContext const &  mc 
)
inherited
bool art::detail::Producer::doEvent ( EventPrincipal &  ep,
ModuleContext const &  mc,
std::atomic< std::size_t > &  counts_run,
std::atomic< std::size_t > &  counts_passed,
std::atomic< std::size_t > &  counts_failed 
)
inherited
void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const &  fb)
inherited
void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const &  fb)
inherited
void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const &  fb)
inherited
void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const &  fb)
inherited
double trk::KalmanTrack::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 
)
private

Definition at line 1755 of file KalmanTrack_module.cc.

References channels, CountHits(), delta, dir, fCellL, fCellW, fDeltaChiAllowed, fNChans, fSliceChans, fSliceChansCellSort, fSysVariance, fView, GetFirstHit(), MECModelEnuComparisons::i, K, PandAna.Demos.pi0_spectra::p0, make_associated_cosmic_defs::p3, elec2geo::pos, R, std::sqrt(), febshutoff_auto::start, ana::weight, and test::z.

Referenced by FindTracks().

1758  {
1759  std::vector<trk::LocatedChan> & channels = zProp? fSliceChans: fSliceChansCellSort;
1760  double chi = 0.0;
1761  double cellVarW = fCellW*fCellW;
1762  double cellVarL = fCellL*fCellL;
1763  double xMinus[3];
1764  double pMinus[3];
1765  double K[2];
1766  int nhits = CountHits(trkHits);
1767  // find the first hit in the track
1768  int firstHit = GetFirstHit(trkHits);
1769  // get the track direction and starting point
1770  const std::array<double, 2> & start = trkTrajs[firstHit];
1771  std::array<double, 2> dir = trkDirs[firstHit];
1772 
1773  // make a vector to store the order hits in (each entry holds a vector of track hits for the given layer)
1774  double trajPointsFor[nhits][3];
1775  double pFor[nhits][4];
1776  double trajPointsBack[nhits][3];
1777  double pBack[nhits][4];
1778  double pos;
1779  double slope;
1780  double otherpos;
1781  // sort all the hits into the order to be added to the track
1782  if(zProp){
1783  pos = start[0];
1784  otherpos = start[1];
1785  if(dir[1] == 0){ dir[1] = 0.0001; }
1786  slope = dir[0]/dir[1];
1787  }
1788  else{
1789  pos = start[1];
1790  otherpos = start[0];
1791  if(dir[0] == 0){ dir[0] = 0.0001; }
1792  slope = dir[1]/dir[0];
1793  }
1794 
1795  std::vector<const rb::CellHit * > trackHitsNew(nhits);
1796  double hitCentersNew[nhits][3];
1797  int ihit = 0;
1798  for(int i = 0; i < fNChans; ++i){
1799  if(trkHits[i]){
1800  trackHitsNew[ihit] = channels[i].GetHit().get();
1801  channels[i].GetCenter(hitCentersNew[ihit]);
1802  ++ihit;
1803  }
1804  }
1805 
1806  // The track is only one plane, no need to filter just choose the cell centers.
1807  if(!zProp && (trackHitsNew[0]->Plane() == trackHitsNew[nhits-1]->Plane())){
1808  // Loop over the trkHits and set trajectories and dirs
1809  int ihit = 0;
1810  for(int i = 0; i < fNChans; ++i){
1811  if(trkHits[i]){
1812  if(zDir > 0){ trkDirs[i][0] = 1.0; }
1813  else{ trkDirs[i][0] = -1.0; }
1814  trkDirs[i][1] = 0.0;
1815  trkTrajs[i][0] = hitCentersNew[ihit][fView];
1816  trkTrajs[i][1] = hitCentersNew[ihit][2];
1817  ++ihit;
1818  }
1819  }
1820  return chi;
1821  }
1822 
1823  // define the measurement and system covariance matrices;
1824  double R = cellVarW/3.0;
1825  if(!zProp){ R = cellVarL/3.0; }
1826  double Q[3];
1827  Q[2] = fSysVariance;
1828 
1829  // filter forward first
1830  int iCurrentHit = 0;
1831  int iNextHit = 0;
1832  int iTraj = 0;
1833  double delta;
1834  // second layer(plane/cell) of hits is at the plane of iNextHit
1835  while(iCurrentHit<nhits && iNextHit<nhits){
1836  int searchPlane = trackHitsNew[iCurrentHit]->Plane();
1837  int searchCell = trackHitsNew[iCurrentHit]->Cell();
1838  int nhitsAdded = 0;
1839  double peTot = 0.0;
1840  if(zProp){
1841  for(int i = iCurrentHit; i<nhits;++i){
1842  if(trackHitsNew[i]->Plane() == searchPlane){
1843  ++nhitsAdded;
1844  ++iNextHit;
1845  peTot+= trackHitsNew[i]->PE();
1846  }
1847  else{ break; }
1848  }
1849  }
1850  else{
1851  for(int i = iCurrentHit; i<nhits; ++i){
1852  if(trackHitsNew[i]->Cell() == searchCell){
1853  ++nhitsAdded;
1854  ++iNextHit;
1855  peTot+= trackHitsNew[i]->PE();
1856  }
1857  else{ break; }
1858  }
1859  }
1860  // loop over all the hits in the layer to find the state estimate of this state.
1861  double otherposOld = otherpos;
1862  double posOld = pos;
1863  double slopeOld = slope;
1864  pos = 0;
1865  otherpos = 0;
1866  slope = 0;
1867  const double pEstimateOld[3] = { pEstimate[0], pEstimate[1], pEstimate[2] };
1868  pEstimate[0] = 0;
1869  pEstimate[1] = 0;
1870  pEstimate[2] = 0;
1871  double z[2];
1872  while(iCurrentHit<iNextHit){
1873  // Get the best hit position from the track hit
1874  if(zProp){
1875  z[0] = hitCentersNew[iCurrentHit][fView];
1876  z[1] = hitCentersNew[iCurrentHit][2];
1877  }
1878  else{
1879  z[0] = hitCentersNew[iCurrentHit][2];
1880  z[1] = hitCentersNew[iCurrentHit][fView];
1881  }
1882  delta = z[1] - otherposOld;
1883 
1884  // update Q and R with new delta
1885  Q[0] = delta*delta*Q[2];
1886  Q[1] = delta*Q[2];
1887 
1888  // Predict next state
1889  xMinus[0] = posOld + delta*slopeOld;
1890  xMinus[2] = slopeOld;
1891  pMinus[0] = pEstimateOld[0] + 2.0*delta*pEstimateOld[1] + delta*delta*pEstimateOld[2] + Q[0];
1892  pMinus[1] = pEstimateOld[1] + delta*pEstimateOld[2] + Q[1];
1893  pMinus[2] = pEstimateOld[2] + Q[2];
1894 
1895  //calculate filter at this layer.
1896  double denom = (pMinus[0]+R);
1897  K[0] = pMinus[0]/denom;
1898  K[1] = pMinus[1]/denom;
1899 
1900  const double weight = trackHitsNew[iCurrentHit]->PE()/peTot;
1901  pos+= weight*(xMinus[0] + K[0]*(z[0]-xMinus[0]));
1902  otherpos+= weight*z[1];
1903  otherpos = z[1];
1904  slope+= weight*(xMinus[2] + K[1]*(z[0]-xMinus[0]));
1905  pEstimate[0]+= weight*(pMinus[0]*R/denom);
1906  pEstimate[1]+= weight*(pMinus[1]*R/denom);
1907  pEstimate[2]+= weight*((pMinus[2]*pMinus[0]-pMinus[1]*pMinus[1]+pMinus[2]*R)/denom);
1908  ++iCurrentHit;
1909  } // end of while loop over iCurrentHit
1910  while(iTraj<iNextHit){
1911  // add the trajectory point and set the new direction at this layer
1912  if(zProp){
1913  trajPointsFor[iTraj][0] = pos;
1914  trajPointsFor[iTraj][1] = otherpos;
1915  }
1916  else{
1917  trajPointsFor[iTraj][0] = otherpos;
1918  trajPointsFor[iTraj][1] = pos;
1919  }
1920  trajPointsFor[iTraj][2] = slope;
1921  pFor[iTraj][0] = pEstimate[0];
1922  pFor[iTraj][3] = pEstimate[2];
1923  ++iTraj;
1924  }
1925  iCurrentHit = iNextHit;
1926  iTraj = iNextHit;
1927  } //end of while loop over hits in the track
1928 
1929  // filter backward
1930  iCurrentHit = nhits-1;
1931  iNextHit = nhits-1;
1932  iTraj = nhits-1;
1933  // second layer(plane/cell) of hits is at the plane of iNextHit
1934  while(iCurrentHit>=0 && iNextHit>=0){
1935  int searchPlane = trackHitsNew[iCurrentHit]->Plane();
1936  int searchCell = trackHitsNew[iCurrentHit]->Cell();
1937  int nhitsAdded = 0;
1938  double peTot = 0.0;
1939  if(zProp){
1940  for(int i = iCurrentHit; i>=0;--i){
1941  if(trackHitsNew[i]->Plane() == searchPlane){
1942  ++nhitsAdded;
1943  --iNextHit;
1944  peTot+= trackHitsNew[i]->PE();
1945  }
1946  else{ break; }
1947  }
1948  }
1949  else{
1950  for(int i = iCurrentHit; i>=0; --i){
1951  if(trackHitsNew[i]->Cell() == searchCell){
1952  ++nhitsAdded;
1953  --iNextHit;
1954  peTot+= trackHitsNew[i]->PE();
1955  }
1956  else{ break; }
1957  }
1958  }
1959  // loop over all the hits in the layer to find the state estimate of this state.
1960  double otherposOld = otherpos;
1961  double posOld = pos;
1962  double slopeOld = slope;
1963  pos = 0;
1964  otherpos = 0;
1965  slope = 0;
1966  const double pEstimateOld[3] = { pEstimate[0], pEstimate[1], pEstimate[2] };
1967  pEstimate[0] = 0;
1968  pEstimate[1] = 0;
1969  pEstimate[2] = 0;
1970  double z[2];
1971  double weight = 1.0/((double)nhitsAdded);
1972  while(iCurrentHit>iNextHit){
1973  // Get the best hit position from the track hit
1974  if(zProp){
1975  z[0] = hitCentersNew[iCurrentHit][fView];
1976  z[1] = hitCentersNew[iCurrentHit][2];
1977  }
1978  else{
1979  z[0] = hitCentersNew[iCurrentHit][2];
1980  z[1] = hitCentersNew[iCurrentHit][fView];
1981  }
1982  delta = z[1] - otherposOld;
1983 
1984  // update Q and R with new delta
1985  Q[0] = delta*delta*Q[2];
1986  Q[1] = delta*Q[2];
1987 
1988  // Predict next state
1989  xMinus[0] = posOld + delta*slopeOld;
1990  xMinus[2] = slopeOld;
1991  pMinus[0] = pEstimateOld[0] + 2.0*delta*pEstimateOld[1] + delta*delta*pEstimateOld[2] + Q[0];
1992  pMinus[1] = pEstimateOld[1] + delta*pEstimateOld[2] + Q[1];
1993  pMinus[2] = pEstimateOld[2] + Q[2];
1994 
1995  //calculate filter at this layer.
1996  double denom = (pMinus[0]+R);
1997  K[0] = pMinus[0]/denom;
1998  K[1] = pMinus[1]/denom;
1999 
2000  weight = trackHitsNew[iCurrentHit]->PE()/peTot;
2001  pos+= weight*(xMinus[0] + K[0]*(z[0]-xMinus[0]));
2002  otherpos+= weight*z[1];
2003  otherpos = z[1];
2004  slope+= weight*(xMinus[2] + K[1]*(z[0]-xMinus[0]));
2005  pEstimate[0]+= weight*(pMinus[0]*R/denom);
2006  pEstimate[1]+= weight*(pMinus[1]*R/denom);
2007  pEstimate[2]+= weight*((pMinus[2]*pMinus[0]-pMinus[1]*pMinus[1]+pMinus[2]*R)/denom);
2008  --iCurrentHit;
2009  } // end of while loop over iCurrentHit
2010  // add the trajectory points
2011  while(iTraj>iNextHit){
2012  if(zProp){
2013  trajPointsBack[iTraj][0] = pos;
2014  trajPointsBack[iTraj][1] = otherpos;
2015  }
2016  else{
2017  trajPointsBack[iTraj][0] = otherpos;
2018  trajPointsBack[iTraj][1] = pos;
2019  }
2020  trajPointsBack[iTraj][2] = slope;
2021  pBack[iTraj][0] = pEstimate[0];
2022  pBack[iTraj][3] = pEstimate[2];
2023  --iTraj;
2024  }
2025  iCurrentHit = iNextHit;
2026  iTraj = iNextHit;
2027  } //end of while loop over hits in the track
2028 
2029  // set trajectory points and calculate the chi^2 of each hit and keep track of the worst outlier
2030  double worstchi = 0;
2031  ihit = 0;
2032  for(int i = 0; i < fNChans; ++i){
2033  if(trkHits[i]){
2034  double trajWeight[2] = {0.5, 0.5};
2035  double dirWeight[2] = {0.5, 0.5};
2036  double p0;
2037  double p3;
2038  double deltaChi = 0.0;
2039  if(zProp){
2040  if(trackHitsNew[ihit]->Plane() == trackHitsNew[0]->Plane()){
2041  trajWeight[0] = 0.0, trajWeight[1] = 1.0;
2042  dirWeight[0] = 0.0, dirWeight[1] = 1.0;
2043  p0 = pBack[ihit][0];
2044  p3 = pBack[ihit][3];
2045  }
2046  else if(trackHitsNew[ihit]->Plane() == trackHitsNew[(nhits-1)]->Plane()){
2047  trajWeight[0] = 1.0, trajWeight[1] = 0.0;
2048  dirWeight[0] = 1.0, dirWeight[1] = 0.0;
2049  p0 = pFor[ihit][0];
2050  p3 = pFor[ihit][3];
2051  }
2052  else{
2053  double trajdenom = 1.0/sqrt(pFor[ihit][0]) + 1.0/sqrt(pBack[ihit][0]);
2054  trajWeight[0] = (1.0/sqrt(pFor[ihit][0]))/trajdenom, trajWeight[1] = (1.0/sqrt(pBack[ihit][0]))/trajdenom;
2055  double dirdenom = 1.0/sqrt(pFor[ihit][3]) + 1.0/sqrt(pBack[ihit][3]);
2056  dirWeight[0] = (1.0/sqrt(pFor[ihit][3]))/dirdenom, dirWeight[1] = (1.0/sqrt(pBack[ihit][3]))/dirdenom;
2057  p0 = 2.0/(trajdenom*trajdenom);
2058  p3 = 2.0/(dirdenom*dirdenom);
2059  }
2060  trkTrajs[i][0] = trajWeight[0]*trajPointsFor[ihit][0] + trajWeight[1]*trajPointsBack[ihit][0];
2061  trkTrajs[i][1] = 0.5*trajPointsFor[ihit][1] + 0.5*trajPointsBack[ihit][1];
2062  slope = dirWeight[0]*trajPointsFor[ihit][2] + dirWeight[1]*trajPointsBack[ihit][2];
2063  double sl2 = slope*slope;
2064  trkDirs[i][0] = slope/sqrt(1.0+sl2);
2065  trkDirs[i][1] = 1.0/sqrt(1.0+sl2);
2066  double r2 = (hitCentersNew[ihit][fView]-trkTrajs[i][0])*(hitCentersNew[ihit][fView]-trkTrajs[i][0]);
2067  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);
2068  }
2069  else{
2070  if(trackHitsNew[ihit]->Cell() == trackHitsNew[0]->Cell()){
2071  trajWeight[0] = 0.0, trajWeight[1] = 1.0;
2072  dirWeight[0] = 0.0, dirWeight[1] = 1.0;
2073  p0 = pBack[ihit][0];
2074  p3 = pBack[ihit][3];
2075  }
2076  else if(trackHitsNew[ihit]->Cell() == trackHitsNew[(nhits-1)]->Cell()){
2077  trajWeight[0] = 1.0, trajWeight[1] = 0.0;
2078  dirWeight[0] = 1.0, dirWeight[1] = 0.0;
2079  p0 = pFor[ihit][0];
2080  p3 = pFor[ihit][3];
2081  }
2082  else{
2083  double trajdenom = 1.0/sqrt(pFor[ihit][0]) + 1.0/sqrt(pBack[ihit][0]);
2084  trajWeight[0] = (1.0/sqrt(pFor[ihit][0]))/trajdenom, trajWeight[1] = (1.0/sqrt(pBack[ihit][0]))/trajdenom;
2085  double dirdenom = 1.0/sqrt(pFor[ihit][3]) + 1.0/sqrt(pBack[ihit][3]);
2086  dirWeight[0] = (1.0/sqrt(pFor[ihit][3]))/dirdenom, dirWeight[1] = (1.0/sqrt(pBack[ihit][3]))/dirdenom;
2087  p0 = 2.0/(trajdenom*trajdenom);
2088  p3 = 2.0/(dirdenom*dirdenom);
2089  }
2090  trkTrajs[i][0] = 0.5*trajPointsFor[ihit][0] + 0.5*trajPointsBack[ihit][0];
2091  trkTrajs[i][1] = trajWeight[0]*trajPointsFor[ihit][1] + trajWeight[1]*trajPointsBack[ihit][1];
2092  slope = dirWeight[0]*trajPointsFor[ihit][2] + dirWeight[1]*trajPointsBack[ihit][2];
2093  double sl2 = slope*slope;
2094  trkDirs[i][0] = 1.0/sqrt(1.0+sl2);
2095  trkDirs[i][1] = slope/sqrt(1.0+sl2);
2096  double r2 = (hitCentersNew[ihit][2]-trkTrajs[i][1])*(hitCentersNew[ihit][2]-trkTrajs[i][1]);
2097  deltaChi+= (r2*(1+sl2)*(1+sl2))/(((1+sl2)*(1+sl2)*(cellVarL/3.0+sl2*cellVarW/3.0+p0))+r2*sl2*p3/2.0);
2098  }
2099  chi+=deltaChi;
2100  // check if this is an outlier
2101  if(deltaChi>fDeltaChiAllowed){
2102  if(deltaChi>worstchi){
2103  worstchi = deltaChi;
2104  newOutlierPos = i;
2105  }
2106  }
2107  chi+=deltaChi;
2108  ++ihit;
2109  } // if hit is on track, ihit
2110  } // end of loop over channels, i
2111 
2112  return chi;
2113  }
delta
double delta
Definition: runWimpSim.h:98
ana::weight
const Var weight
Definition: PIDCutTuningDefs.h:96
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
make_associated_cosmic_defs.p3
p3
Definition: make_associated_cosmic_defs.py:192
K
Double_t K
Definition: TG43_relative_dose.C:17
febshutoff_auto.start
start
Definition: febshutoff_auto.py:53
trk::KalmanTrack::fSliceChansCellSort
std::vector< trk::LocatedChan > fSliceChansCellSort
Definition: KalmanTrack_module.cc:66
channels
std::map< ToFCounter, std::vector< unsigned int > > channels
Definition: ToFPositionRecoAnalysis_module.cc:81
trk::KalmanTrack::fSliceChans
std::vector< trk::LocatedChan > fSliceChans
Definition: KalmanTrack_module.cc:65
R
#define R(x)
test.z
z
Definition: test.py:28
PandAna.Demos.pi0_spectra.p0
p0
Definition: pi0_spectra.py:151
dir
TDirectory * dir
Definition: macro.C:5
trk::KalmanTrack::CountHits
int CountHits(const std::vector< bool > &trkHits)
Definition: KalmanTrack_module.cc:2468
trk::KalmanTrack::GetFirstHit
int GetFirstHit(const std::vector< bool > &trkHits)
Definition: KalmanTrack_module.cc:2503
double trk::KalmanTrack::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 
)
private

Definition at line 1104 of file KalmanTrack_module.cc.

References abs(), trk::KalmanGeoHelper::BestTrackPoint(), channels, CountHits(), trk::KalmanGeoHelper::CountMissedCellsOnLine(), delta, geo::GeometryBase::DetHalfHeight(), geo::GeometryBase::DetHalfWidth(), geo::GeometryBase::DetLength(), dir, fbadc, fBadChannels, fCellL, fCellW, fDeltaChiAllowed, fGapAllowed, fgeom, fNChans, fSliceChans, fSliceChansCellSort, fSysVariance, fView, MECModelEnuComparisons::i, chaninfo::BadChanList::IsBad(), chaninfo::BadChanList::IsLowEfficiency(), K, geo::kPLANE_NOT_FOUND, NCells, geo::GeometryBase::NextPlaneInView(), PandAna.reco_validation.add_data::offset, makeBrightnessMap::position, R, std::sqrt(), febshutoff_auto::start, makeDatasetsPage::state, and plot_timing_data::swap.

Referenced by FindTracks().

1110  {
1111  std::vector<trk::LocatedChan> & channels = zProp? fSliceChans : fSliceChansCellSort;
1112  double chi = 0.0; // return value
1113  // get which 2d view we are in from the slice hits
1114  double maxDeltaChi = fDeltaChiAllowed;
1115  int maxMisses = 1.5*fGapAllowed;
1116  double cellVarW = fCellW*fCellW;
1117  double cellVarL = fCellL*fCellL;
1118  int misses = 0;
1119  bool done = false;
1120  int nHits = CountHits(trkHits);
1121  if(nHits == 0){ return 0.0; }
1122 
1123  // find the first hit attached to this track
1124  int iFirstHit = 0;
1125 
1126  if(zDir < 0){ iFirstHit = fNChans-1; }
1127  if(firsthit < 0){
1128  bool firstHitfound = false;
1129  while(!firstHitfound && iFirstHit<fNChans && iFirstHit >= 0){
1130  if(zProp){
1131  if(trkHits[iFirstHit]){ firstHitfound = true; }
1132  else{iFirstHit+=zDir;}
1133  }
1134  else{
1135  if(trkHits[iFirstHit]){ firstHitfound = true; }
1136  else{ iFirstHit+=zDir; }
1137  }
1138  }
1139  }
1140  else{ iFirstHit = firsthit; }
1141 
1142  unsigned short iPlane = channels[iFirstHit].Plane();
1143  unsigned short iCell = channels[iFirstHit].Cell();;
1144 
1145  // Find the first slice hit that belongs to the same plane/cell as the first hit attached to the track
1146  int iCurrentHit = 0;
1147  if(zDir < 0){ iCurrentHit = fNChans-1; }
1148  bool currentHitfound = false;
1149  while(!currentHitfound && iCurrentHit<fNChans && iCurrentHit >= 0){
1150  if(zProp){
1151  if(channels[iCurrentHit].Plane() == iPlane){ currentHitfound = true; }
1152  else{iCurrentHit+=zDir;}
1153  }
1154  else{
1155  if(channels[iCurrentHit].Cell() == iCell){ currentHitfound = true; }
1156  else{ iCurrentHit+=zDir; }
1157  }
1158  }
1159 
1160  // get the starting trajectory point of the track (trajectory point of the first hit)
1161  const std::array<double, 2> & start = trkTrajs[iFirstHit];
1162  std::array<double, 2> dir = trkDirs[iFirstHit];
1163  int iCurrentTrackHit = iFirstHit; // keep track of which index corresponds to the currentHitPlane, currentHitCell
1164 
1165  double position;
1166  double slope;
1167  double otherpos;
1168  if(zProp){
1169  position = start[0];
1170  otherpos = start[1];
1171  if(dir[1] == 0){ dir[1] = 0.0001; }
1172  slope = dir[0]/dir[1];
1173  }
1174  else{
1175  position = start[1];
1176  otherpos = start[0];
1177  if(dir[0] == 0){ dir[0] = 0.0001; }
1178  slope = dir[1]/dir[0];
1179  }
1180 
1181  double R = cellVarW/3.0;
1182  if(!zProp){ R = cellVarL/3.0; };
1183  double Q[3] = {0, 0, fSysVariance};
1184  // state[0] = filtered position, state[1] = position of measurement in other view, state[2] = slope;
1185  double state[3] = {position, otherpos, slope};
1186 
1187  int iNextHit = iCurrentHit;
1188 
1189  // forget about all current track hits
1190  if(firsthit < 0){
1191  for(int i = 0; i < fNChans; ++i){
1192  trkHits[i] = false;
1193  }
1194  outnhits = 0;
1195  }
1196  else{ outnhits = nHits; }
1197 
1198  int nTrackHits = outnhits; // running total of number of hits added to track
1199  // Create variable for projecting states into the next plane of the detector
1200  // initialize the estimates state and error covairance
1201  double xMinus[3];
1202  double pMinus[3];
1203  double K[2];
1204  double xPredict[3];
1205  double cPredict[3];
1206 
1207  double delta;
1208  short currentHitPlane = iPlane;
1209  int expectedHitPlane = currentHitPlane;
1210  short currentHitCell = iCell;
1211  short expectedHitCell = currentHitCell;
1212 
1213  std::array<double, 2> bestpos;
1214 
1215  while(((iCurrentHit < (fNChans-1) && zDir > 0) || (iCurrentHit > 0 && zDir < 0)) && !done){
1216 
1217  short nextHitPlane = channels[iNextHit].Plane();
1218  short nextHitCell = channels[iNextHit].Cell();
1219 
1220  // update current hit plane and cell
1221  if(nTrackHits>0){
1222  currentHitPlane = channels[iCurrentTrackHit].Plane();
1223  currentHitCell = channels[iCurrentTrackHit].Cell();
1224  }
1225  // Define a boolean to check if the next hit is in the layer(plane/cell) that we expect it to be in;
1226  bool sameLayerTest;
1227  if(zProp){sameLayerTest = ((nextHitPlane-currentHitPlane) == (expectedHitPlane-currentHitPlane)); }
1228  else{ sameLayerTest = (((nextHitCell-currentHitCell) == (expectedHitCell-currentHitCell))); }
1229 
1230  // Check if the next slice hit is in the plane/cell where the expected hit
1231  if(sameLayerTest){
1232  // keep track of information about hits added to the track
1233  double peWeight = 0;
1234  double layerState[3] = {state[0], state[1], state[2]};
1235  double layerP[3] = {errorCov[0], errorCov[1], errorCov[2]};
1236 
1237  int k = iNextHit;
1238  int firstLayerHit = iNextHit;
1239 
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;
1250  --outnhits;
1251  }
1252  iCurrentHit = iNextHit;
1253  iNextHit = iCurrentHit + zDir;
1254  }
1255  bool doneAdding = false;
1256  if(nadd == 0){ doneAdding = true;}
1257 
1258  int nhits = 0;
1259  int layerWindowlow = 0;
1260  int layerWindowhigh = 999999;
1261  if(!zProp){
1262  if(((state[2] <= 0 && zDir < 0) || (state[2] >= 0 && zDir > 0))){
1263  if(nTrackHits>0){
1264  layerWindowlow = channels[iCurrentTrackHit].Plane();
1265  if(abs(nextHitCell-currentHitCell)>1){
1266  layerWindowlow = channels[iCurrentTrackHit].Plane()-2;
1267  layerWindowhigh = channels[iCurrentTrackHit].Plane()+6;
1268  }
1269  }
1270  }
1271  else{
1272  if(nTrackHits>0){
1273  layerWindowhigh = channels[iCurrentTrackHit].Plane();
1274  if(abs(nextHitCell-currentHitCell)>1){
1275  layerWindowhigh = channels[iCurrentTrackHit].Plane()+2;
1276  layerWindowlow = channels[iCurrentTrackHit].Plane()-6;
1277  }
1278  }
1279  }
1280  }
1281  while(!doneAdding){
1282  double minDistBest = -1;
1283  int bestHitIndex = 9999;
1284  if(nadd == 1 && !layerTrackHits[0]){ bestHitIndex = 0; }
1285  else{
1286  for(int i = 0; i<nadd;++i){
1287  // do nothing if this hit has already been added to the track
1288  if(!layerTrackHits[i]){
1289  int currChanIdx = firstLayerHit+i*zDir;
1290 
1291  // check if this hit falls within the track window
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; }
1294 
1295  //Get the bestpoint of this track
1296  double xyzHit[3];
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]);}
1300  // Get the minimum distance to the projected track
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;
1307  }
1308  }
1309  } // end of loop over nadd
1310  } // end of else nadd == 1
1311  if(bestHitIndex != 9999){
1312  int bestChanIdx = firstLayerHit+bestHitIndex;
1313  double xyzHit[3];
1314  channels[bestChanIdx].GetCenter(xyzHit);
1315 
1316  // Update the prediction and error matrices
1317  if(zProp){ delta = xyzHit[2] - layerState[1]; }
1318  else{ delta = xyzHit[fView] - layerState[1]; }
1319  Q[0] = delta*delta*Q[2] + layerP[0];
1320  Q[1] = delta*Q[2];
1321 
1322  //Predict next state using a linear model
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];
1328 
1329  //calculate filter at this layer
1330  double denom = (pMinus[0]+R);
1331  K[0] = pMinus[0]/denom;
1332  K[1] = pMinus[1]/denom;
1333 
1334  // get the measurements from the hit
1335  double candidatePos;
1336  double ipos;
1337  if(zProp){
1338  candidatePos = xyzHit[fView];
1339  ipos = xyzHit[2];
1340  xPredict[1] = ipos;
1341  }
1342  else{
1343  candidatePos = xyzHit[2];
1344  ipos = xyzHit[fView];
1345  xPredict[1] = ipos;
1346  }
1347 
1348  //calculate corrected states with inclusion of this hit
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;
1354 
1355  double r2 = (candidatePos - xPredict[0])*(candidatePos - xPredict[0]);
1356  double sl2 = xPredict[2]*xPredict[2];
1357  double deltaChi = 0.0;
1358  if(zProp){
1359  deltaChi = (r2*(1+sl2)*(1+sl2))/
1360  (((1+sl2)*(1+sl2)*(cellVarW/3+sl2*cellVarL/3+cPredict[0]))+r2*sl2*cPredict[2]/2);
1361  }
1362  else{
1363  deltaChi = (r2*(1+sl2)*(1+sl2))/
1364  (((1+sl2)*(1+sl2)*(cellVarL/3+sl2*cellVarW/3+cPredict[0]))+r2*sl2*cPredict[2]/2);
1365  }
1366 
1367  if(deltaChi <= maxDeltaChi){
1368  ++nhits;
1369  chi+=deltaChi;
1370  trkHits[bestChanIdx] = true;
1371  ++outnhits;
1372  ++nTrackHits;
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;
1378  if(nhits == 1){
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;
1386  }
1387  else{
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;
1395  }
1396  if(zProp){
1397  // update dir
1398  dir[0] = layerState[2]/sqrt(1.0+layerState[2]*layerState[2]);
1399  dir[1] = 1.0/sqrt(1.0+layerState[2]*layerState[2]);
1400  if(nhits == 1){
1401  int newLayerWindowlow = currentHitCell-fGapAllowed;
1402  int newLayerWindowhigh = currentHitCell+fGapAllowed;
1403  if(fBadChannels){
1404  while(newLayerWindowlow > 0 &&
1405  (fbadc->IsBad(currentHitPlane,newLayerWindowlow)
1406  || fbadc->IsLowEfficiency(currentHitPlane,newLayerWindowlow))){
1407  --newLayerWindowlow;
1408  }
1409  while(newLayerWindowhigh < NCells[fView] &&
1410  (fbadc->IsBad(currentHitPlane,newLayerWindowhigh)
1411  || fbadc->IsLowEfficiency(currentHitPlane,newLayerWindowhigh))){
1412  ++newLayerWindowhigh;
1413  }
1414  }
1415  layerWindowlow = newLayerWindowlow;
1416  layerWindowhigh = newLayerWindowhigh;
1417  }
1418  else{
1419  if((currentHitCell - layerWindowlow) < (layerWindowhigh - currentHitCell)){
1420  if(!fBadChannels){layerWindowlow = currentHitCell-fGapAllowed;}
1421  else{
1422  // increase layer window by 1 cell until no more bad channels found
1423  int newLayerWindowlow = currentHitCell-fGapAllowed;
1424  while(newLayerWindowlow > 0 &&
1425  (fbadc->IsBad(currentHitPlane,newLayerWindowlow)
1426  || fbadc->IsLowEfficiency(currentHitPlane,newLayerWindowlow)) ){
1427  --newLayerWindowlow;
1428  }
1429  layerWindowlow = newLayerWindowlow;
1430  }
1431  }
1432  else{
1433  if(!fBadChannels){layerWindowhigh = currentHitCell+fGapAllowed;}
1434  else{
1435  // increase layer window by 1 cell until no more bad channels found
1436  int newLayerWindowhigh = currentHitCell+fGapAllowed;
1437  while(newLayerWindowhigh < NCells[fView] &&
1438  (fbadc->IsBad(currentHitPlane,newLayerWindowhigh)
1439  || fbadc->IsLowEfficiency(currentHitPlane,newLayerWindowhigh)) ){
1440  ++newLayerWindowhigh;
1441  }
1442  layerWindowhigh = newLayerWindowhigh;
1443  }
1444  }
1445  }
1446  } // end of if zProp
1447  else{
1448  // update dir
1449  dir[0] = 1.0/sqrt(1.0+layerState[2]*layerState[2]);
1450  dir[1] = layerState[2]/sqrt(1.0+layerState[2]*layerState[2]);
1451  if(nhits == 1){
1452  if(((state[2] <= 0 && zDir < 0) || (state[2] >= 0 && zDir > 0))){
1453  int newLayerWindowlow = currentHitPlane;
1454  int newLayerWindowhigh = fgeom->NextPlaneInView(currentHitPlane,+1);
1455  layerWindowlow = newLayerWindowlow;
1456  if(newLayerWindowhigh == geo::kPLANE_NOT_FOUND){newLayerWindowhigh = currentHitPlane;}
1457  else{
1458  if(fBadChannels){
1459  int newlayerWindowhightest = newLayerWindowhigh;
1460  while(newlayerWindowhightest != geo::kPLANE_NOT_FOUND &&
1461  (fbadc->IsBad(newlayerWindowhightest,currentHitCell)
1462  || fbadc->IsLowEfficiency(newlayerWindowhightest,currentHitCell)) ){
1463  newlayerWindowhightest = fgeom->NextPlaneInView(newLayerWindowhigh,+1);
1464  if(newlayerWindowhightest != geo::kPLANE_NOT_FOUND){newLayerWindowhigh = newlayerWindowhightest;}
1465  }
1466  }
1467  }
1468  layerWindowhigh = newLayerWindowhigh;
1469  }
1470  else{
1471  int newLayerWindowlow = fgeom->NextPlaneInView(currentHitPlane,-1);
1472  int newLayerWindowhigh = currentHitPlane;
1473  layerWindowhigh = newLayerWindowhigh;
1474  if(newLayerWindowlow == geo::kPLANE_NOT_FOUND){newLayerWindowlow = currentHitPlane;}
1475  else{
1476  if(fBadChannels){
1477  int newlayerWindowlowtest = newLayerWindowlow;
1478  while(newlayerWindowlowtest != geo::kPLANE_NOT_FOUND &&
1479  (fbadc->IsBad(newlayerWindowlowtest,currentHitCell)
1480  || fbadc->IsLowEfficiency(newlayerWindowlowtest,currentHitCell)) ){
1481  newlayerWindowlowtest = fgeom->NextPlaneInView(newLayerWindowlow,-1);
1482  if(newlayerWindowlowtest != geo::kPLANE_NOT_FOUND){newLayerWindowlow = newlayerWindowlowtest;}
1483  }
1484  }
1485  }
1486  layerWindowlow = newLayerWindowlow;
1487  }
1488  }
1489  else{
1490  if(((state[2] <= 0 && zDir < 0) || (state[2] >= 0 && zDir > 0))){
1491  layerWindowlow = currentHitPlane;
1492  if(!fBadChannels){
1493  layerWindowhigh = fgeom->NextPlaneInView(currentHitPlane,1);
1494  if(layerWindowhigh == geo::kPLANE_NOT_FOUND){layerWindowhigh = currentHitPlane;}
1495  }
1496  else{
1497  int newLayerWindowhigh = fgeom->NextPlaneInView(currentHitPlane,1);
1498  if(newLayerWindowhigh == geo::kPLANE_NOT_FOUND){newLayerWindowhigh = currentHitPlane;}
1499  else{
1500  int newLayerWindowhightest = newLayerWindowhigh;
1501  while(newLayerWindowhightest != geo::kPLANE_NOT_FOUND &&
1502  (fbadc->IsBad(newLayerWindowhightest,currentHitCell)
1503  || fbadc->IsLowEfficiency(newLayerWindowhightest,currentHitCell)) ){
1504  newLayerWindowhightest = fgeom->NextPlaneInView(newLayerWindowhigh,+1);
1505  if(newLayerWindowhightest != geo::kPLANE_NOT_FOUND){newLayerWindowhigh = newLayerWindowhightest;}
1506  }
1507  }
1508  layerWindowhigh = newLayerWindowhigh;
1509  }
1510  }
1511  else{
1512  layerWindowhigh = currentHitPlane;
1513  if(!fBadChannels){
1514  layerWindowlow = fgeom->NextPlaneInView(currentHitPlane,-1);
1515  if(layerWindowlow == geo::kPLANE_NOT_FOUND){layerWindowlow = currentHitPlane;}
1516  }
1517  else{
1518  int newLayerWindowlow = fgeom->NextPlaneInView(currentHitPlane,-1);
1519  if(newLayerWindowlow == geo::kPLANE_NOT_FOUND){newLayerWindowlow = currentHitPlane;}
1520  else{
1521  int newLayerWindowlowtest = newLayerWindowlow;
1522  while(newLayerWindowlowtest != geo::kPLANE_NOT_FOUND &&
1523  (fbadc->IsBad(newLayerWindowlowtest,currentHitCell)
1524  || fbadc->IsLowEfficiency(newLayerWindowlowtest,currentHitCell)) ){
1525  newLayerWindowlowtest = fgeom->NextPlaneInView(newLayerWindowlow,-1);
1526  if(newLayerWindowlowtest != geo::kPLANE_NOT_FOUND){newLayerWindowlow = newLayerWindowlowtest;}
1527  }
1528  }
1529  layerWindowlow = newLayerWindowlow;
1530  }
1531  }
1532  }
1533  } // end of else of if(zProp) else
1534  } // end of if(deltaChi <= maxDeltaChi)
1535  else{ layerTrackHits[(bestHitIndex/zDir)] = true; }
1536  // now we know everything about this hit delete it from our knowledge bank
1537  } // if(bestHitIndex != 9999)
1538  else{ doneAdding = true;}
1539  } // while(!doneAdding)
1540 
1541  // update the filtered state and covariance matrix based on added hits
1542  if(nhits>0){
1543  // Add track information to output variables
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];
1550 
1551  // update track and find next expected hit layer
1552  if(zProp){
1553  if(fLayerHit < iLayerHit){ std::swap(iLayerHit,fLayerHit); }
1554  for(int ihit = iLayerHit; ihit<=fLayerHit; ++ihit){
1555  if(trkHits[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];
1560  }
1561  }
1562  expectedHitPlane = fgeom->NextPlaneInView(currentHitPlane,zDir);
1563  if(expectedHitPlane == geo::kPLANE_NOT_FOUND){ done = true; }
1564  }
1565  else{
1566  if(fLayerHit < iLayerHit){ std::swap(iLayerHit,fLayerHit); }
1567  for(int ihit = iLayerHit; ihit<=fLayerHit; ++ihit){
1568  if(trkHits[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];
1573  }
1574  }
1575  expectedHitCell = currentHitCell+zDir;
1576  if(expectedHitCell < 0 || expectedHitCell > NCells[fView] ){ done = true; }
1577  }
1578  if(iNextHit > (fNChans-1) || iNextHit < 0){ done = true; }
1579  // reset number of misses to zero
1580  misses = 0;
1581  } // if(nhits > 0)
1582  else{
1583  // check if there could possibly be anymore hits to add to the track
1584  if(iNextHit > (fNChans-1) || iNextHit < 0){ done = true; }
1585  else{
1586  short nextHitPlane = channels[iNextHit].Plane();
1587  short nextHitCell = channels[iNextHit].Cell();
1588  // project the last trajectory point to the beginning of the expected hit plane
1589  double currtraj[2]; // [0] = v, [1] = z
1590  //point on next cell that track will hit
1591  double nextpt[3]; // [0] = x, [1] = y, [2] = z
1592  channels[iNextHit].GetCenter(nextpt);
1593  if(zProp){
1594  // change current trajectory to just outside of the current hit plane
1595  double offset = zDir*(fCellL+2.0);
1596  currtraj[1] = state[1] + offset;
1597  currtraj[0] = state[2]*offset + state[0];
1598  // change next z position to the z location just outside of the next hit plane
1599  nextpt[2]-=offset;
1600  nextpt[fView] = state[2]*(nextpt[2]-currtraj[1])+currtraj[0];
1601  // check if this point is out of the detector
1602  if(fView == 0 && abs(nextpt[fView]) > fgeom->DetHalfWidth()){ done = true; }
1603  else if(fView == 1 && abs(nextpt[fView]) > fgeom->DetHalfHeight()){ done = true; }
1604  }
1605  else{
1606  // change current trajectory to just outside of the current hit cell
1607  double offset = zDir*(fCellW+2.0);
1608  currtraj[0] = state[1] + offset;
1609  currtraj[1] = state[2]*offset + state[0];
1610  // change next hit position to just outside of the next hit cell
1611  nextpt[fView]-=offset;
1612  nextpt[2] = state[2]*(nextpt[fView]-currtraj[0])+currtraj[1];
1613  if(nextpt[2] < 0 || nextpt[2] > fgeom->DetLength()){ done = true; }
1614  }
1615  misses+=kgeo.CountMissedCellsOnLine(fView,currtraj[0],currtraj[1],nextpt[fView],nextpt[2],fGapAllowed);
1616  if(misses > maxMisses){done = true;}
1617  else{
1618  expectedHitCell = nextHitCell;
1619  expectedHitPlane = nextHitPlane;
1620  }
1621  }
1622  } // end of else for checking if more hits exist in slice past this one
1623  } // end if in samelayertest
1624  else{
1625  // check how many cells this track goes through before getting to the next hit plane
1626  double currtraj[2]; // [0] = v, [1] = z
1627  //point on next cell that track will hit
1628  double nextpt[3]; // [0] = x, [1] = y, [2] = z
1629  channels[iNextHit].GetCenter(nextpt);
1630  if(zProp){
1631  // change current trajectory to just outside of the current hit plane
1632  double offset = zDir*(fCellL+2.0);
1633  currtraj[1] = state[1] + offset;
1634  currtraj[0] = state[2]*offset + state[0];
1635  // change next z position to the z location just outside of the next hit plane
1636  nextpt[2]-=offset;
1637  nextpt[fView] = state[2]*(nextpt[2]-currtraj[1])+currtraj[0];
1638  // check if this point is out of the detector
1639  if(fView == 0 && abs(nextpt[fView]) > fgeom->DetHalfWidth()){ done = true; }
1640  else if(fView == 1 && abs(nextpt[fView]) > fgeom->DetHalfHeight()){ done = true; }
1641  }
1642  else{
1643  // change current trajectory to just outside of the current hit cell
1644  double offset = zDir*(fCellW+2.0);
1645  currtraj[0] = state[1] + offset;
1646  currtraj[1] = state[2]*offset + state[0];
1647  // change next hit position to just outside of the next hit cell
1648  nextpt[fView]-=offset;
1649  nextpt[2] = state[2]*(nextpt[fView]-currtraj[0])+currtraj[1];
1650  if(nextpt[2] < 0 || nextpt[2] > fgeom->DetLength()){ done = true; }
1651  }
1652  misses+=kgeo.CountMissedCellsOnLine(fView,currtraj[0],currtraj[1],nextpt[fView],nextpt[2],fGapAllowed);
1653  if(misses > maxMisses){done = true;}
1654  else{
1655  expectedHitPlane = nextHitPlane;
1656  expectedHitCell = nextHitCell;
1657  }
1658 
1659  } //end of else in samelayertest
1660 
1661  } // end of while loop over all slicehits
1662 
1663  return chi;
1664  }
delta
double delta
Definition: runWimpSim.h:98
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
geo::GeometryBase::DetLength
double DetLength() const
Definition: GeometryBase.cxx:481
abs
void abs(TH1 *hist)
Definition: wrong_sign_uncertainty.C:716
K
Double_t K
Definition: TG43_relative_dose.C:17
febshutoff_auto.start
start
Definition: febshutoff_auto.py:53
trk::KalmanTrack::fSliceChansCellSort
std::vector< trk::LocatedChan > fSliceChansCellSort
Definition: KalmanTrack_module.cc:66
channels
std::map< ToFCounter, std::vector< unsigned int > > channels
Definition: ToFPositionRecoAnalysis_module.cc:81
plot_timing_data.swap
string swap
Definition: plot_timing_data.py:51
trk::KalmanTrack::fbadc
art::ServiceHandle< chaninfo::BadChanList > fbadc
Definition: KalmanTrack_module.cc:55
trk::KalmanTrack::fSliceChans
std::vector< trk::LocatedChan > fSliceChans
Definition: KalmanTrack_module.cc:65
R
#define R(x)
geo::GeometryBase::DetHalfHeight
double DetHalfHeight() const
Definition: GeometryBase.cxx:474
PandAna.reco_validation.add_data.offset
offset
Definition: add_data.py:36
geo::GeometryBase::DetHalfWidth
double DetHalfWidth() const
Definition: GeometryBase.cxx:467
dir
TDirectory * dir
Definition: macro.C:5
trk::KalmanTrack::fgeom
art::ServiceHandle< geo::Geometry > fgeom
Definition: KalmanTrack_module.cc:54
geo::GeometryBase::NextPlaneInView
const unsigned int NextPlaneInView(unsigned int p, int d=+1) const
Definition: GeometryBase.cxx:936
trk::KalmanTrack::CountHits
int CountHits(const std::vector< bool > &trkHits)
Definition: KalmanTrack_module.cc:2468
chaninfo::BadChanList::IsBad
bool IsBad(int plane, int cell)
Definition: BadChanList_service.cc:563
chaninfo::BadChanList::IsLowEfficiency
bool IsLowEfficiency(int plane, int cell)
Definition: BadChanList_service.cc:691
std::vector< rb::Track > trk::KalmanTrack::FindTracks ( KalmanGeoHelper &  kgeo,
std::vector< trk::LocatedChan >  sliceChans 
)
private

Definition at line 399 of file KalmanTrack_module.cc.

References abs(), rb::Cluster::Add(), rb::Track::AppendTrajectoryPoint(), avgX, avgY, trk::KalmanGeoHelper::BestTrackPoint(), rb::CellHit::Cell(), rb::Cluster::Cell(), CellToPlaneSort(), CellToPlaneSortHits(), CheckTrack(), rb::Track::ClearTrajectoryPoints(), CountHits(), CreateHitMaps(), dir, rb::Prong::Dir(), fgeom, FilterOnly(), FilterTracker(), fLargeSliceHits, fLongTrackLoop, fMinHits, fNChans, fSliceChans, fSliceChansCellSort, fTrySingle, fView, geo::CellGeo::GetCenter(), GetFirstHit(), geo::CellGeo::HalfD(), geo::CellGeo::HalfW(), MECModelEnuComparisons::i, IsSinglePlane(), it, rb::Cluster::MaxPlane(), rb::Cluster::MinPlane(), rb::Cluster::NCell(), rb::Track::NTrajectoryPoints(), rb::CellHit::Plane(), geo::GeometryBase::Plane(), PlaneToCellSort(), PlaneToCellSortHits(), RemoveHitsFromSignal(), RemoveSameTracks(), rb::Track::RemoveTrajectoryPoint(), site_stats_from_log::reverse, SegmentFinder(), rb::Track::SetStart(), trk::SortByCell(), trk::SortByPlane(), febshutoff_auto::start, rb::Prong::Start(), rb::Track::Stop(), rb::Track::TrajectoryPoint(), X, Y, and Z.

Referenced by produce().

400  {
401  std::vector<rb::Track> outtracks;
402  // set a boolean to tell when we are done finding tracks
403  bool Done = false;
404  // keep track of times in the loop so only the first time look for long tracks only - SL
405  int nTimes = 0;
406  fTrySingle = true;
407  // Loop over the slice hits until no more acceptable tracks are found
408  // Now that we have cell hits lets sort them from lowest plane to highest
409 
410  trk::SortByPlane(fSliceChans);
411  while(!Done){
412  fSliceChansCellSort = fSliceChans;
413  trk::SortByCell(fSliceChansCellSort);
414  // (re)create hit map between plane and cell sorted channels
415  CreateHitMaps();
416 
417  // Make segments which consist of two points that form the start of a track
418  std::vector<rb::Track> trackSegs;
419  std::vector< std::vector<double> > iP;
420  std::vector<bool> zProps;
421  std::vector< std::vector<bool> > trkHits;
422  std::vector< std::vector<std::array<double, 2> > > trkTrajs;
423  std::vector< std::vector<std::array<double, 2> > > trkDirs;
424  bool initTrack = true;
425  if(fNChans < fLargeSliceHits){ fTrySingle = false; }
426 
427  int ntracks = SegmentFinder(kgeo,iP,zProps,trkHits,trkTrajs,trkDirs);
428 
429  //create some vectors to hold track information in
430  std::vector< double > Chi2(trkHits.size(),0.0); // chi squared value of the track;
431  std::vector< int > nTrkHits(trkHits.size()); // number of hits in the track
432  std::vector< bool > ignoreTrk(trkHits.size(),false); // ignore this track?
433 
434  unsigned int propDirIters = 3; // maximum number of iterations to do tracking/filtering
435  int propDir = -1;
436  for(unsigned int iprop = 0; iprop < propDirIters; ++iprop){
437 
438  // Do the track finding
439  int largestTrack = 0;
440  for(unsigned int it = 0; it< trkHits.size(); ++it){
441  if(!ignoreTrk[it]){
442  int n = CountHits(trkHits[it]);
443  nTrkHits[it] = n;
444  if(n > 0 && n < fNChans){
445  if(n > 2 && fTrySingle){
446  Chi2[it] = FilterTracker(kgeo,propDir,zProps[it],trkHits[it],trkTrajs[it],trkDirs[it],iP[it],nTrkHits[it]);
447  }
448  else{
449  Chi2[it] = FilterTracker(kgeo,propDir,zProps[it],trkHits[it],trkTrajs[it],trkDirs[it],iP[it],nTrkHits[it]);
450  }
451  if(nTrkHits[it] > largestTrack){ largestTrack = nTrkHits[it]; }
452  }
453  }
454  }
455  if(fTrySingle && ntracks > 0 && nTrkHits[0] == fNChans){ initTrack = false; }
456 
457  // Remove tracks that have exactly the same hits in them as another track, keeping the one with the lowest chi square.
458  RemoveSameTracks(trkHits,trkTrajs,trkDirs,Chi2,iP,zProps,nTrkHits,ignoreTrk);
459 
460  propDir = (-1*propDir); // Set filter direction to opposite of that in the track finding stage
461  int minFiltHits = 2; // Need at least two hits to make a track
462 
463  // If we're on the last propogation stage, then minimum filtered hits
464  // is the minimum number of hits for a track
465  if(iprop == propDirIters-1 && fMinHits>2){ minFiltHits = (fMinHits-1); }
466 
467  // Iterate over the filtering
468  for(unsigned int i = 0; i< trkHits.size(); i++){
469 
470  // no need to filter tracks that we don't care about
471  if(nTrkHits[i] >= minFiltHits && !ignoreTrk[i]){
472 
473  // check that this is being propogated in the correct manner
474  if(zProps[i]){
475  if(IsSinglePlane(trkHits[i])){
476  zProps[i] = false;
477  PlaneToCellSortHits(trkHits[i]);
478  PlaneToCellSort(trkTrajs[i]);
479  PlaneToCellSort(trkDirs[i]);
480  }
481  }
482  else{
483  std::vector<bool> sortedhits = trkHits[i];
484  CellToPlaneSortHits(sortedhits);
485  if(!IsSinglePlane(sortedhits)){
486  zProps[i] = true;
487  trkHits[i] = sortedhits;
488  CellToPlaneSort(trkTrajs[i]);
489  CellToPlaneSort(trkDirs[i]);
490  }
491  }
492 
493  // define a boolean to determine whether or not to continue iterating
494  bool filter = true;
495  // save the number of hits on the track
496  int nhitsold = nTrkHits[i];
497 
498  // keep track of the number of iterations needed for filtering
499  int filterIter = 0;
500  int maxFilter = 2; // maximum number of iterations over filtering
501 
502  while(filter && filterIter < maxFilter){
503 
504  ++filterIter; // update the counter
505 
506  // Filter the track
507  int trkOutlier = -1;
508  Chi2[i] = FilterOnly(trkHits[i],trkTrajs[i],trkDirs[i],iP[i],propDir,zProps[i],trkOutlier);
509 
510  // remove outlier if it exists, give the fitter at least one chance to get it right first
511  if(filterIter>0){
512  if(trkOutlier != -1){
513  trkHits[i][trkOutlier] = false;
514  --nTrkHits[i];
515  --filterIter;
516  }
517  }
518 
519  // if number of hits is less than the minimum than stop filtering
520  if(nTrkHits[i] < minFiltHits){
521  filter = false;
522  break; // exit while loop of filtering
523  }
524 
525  // May need to re-sort hits if any where removed
526  if(nTrkHits[i] != nhitsold){
527  if(zProps[i]){
528  if(IsSinglePlane(trkHits[i])){
529  zProps[i] = false;
530  PlaneToCellSortHits(trkHits[i]);
531  PlaneToCellSort(trkTrajs[i]);
532  PlaneToCellSort(trkDirs[i]);
533  }
534  }
535  }
536 
537  // If we didn't remove any hits and filtered more than once we can stop
538  if(nhitsold == nTrkHits[i] && filterIter > 0){
539 
540  // If this is the last propogation
541  // check that this track meets the tracking criteria
542  if(iprop == propDirIters-1 || !initTrack){
543 
544  // check the track
545  std::vector<bool> newTrkHits;
546  CheckTrack(kgeo,trkHits[i],trkTrajs[i],zProps[i],newTrkHits);
547 
548  // update hit counters
549  nTrkHits[i] = CountHits(trkHits[i]);
550  int nhitsnew = CountHits(newTrkHits);
551 
552  if(nhitsnew > 0 && nTrkHits[i] > 0){
553  --filterIter;
554 
555  // if track now has less than desired number of hits, stop filtering
556  if(nTrkHits[i] < minFiltHits){ filter = false; }
557 
558  // may need to re-sort hits if this track had hits removed
559  if(filter && zProps[i]){
560  if(IsSinglePlane(trkHits[i])){
561  zProps[i] = false;
562  PlaneToCellSortHits(trkHits[i]);
563  PlaneToCellSort(trkTrajs[i]);
564  PlaneToCellSort(trkDirs[i]);
565  }
566  }
567 
568  // add the second piece of the track onto the track list if its large enough
569  if(nhitsnew >= minFiltHits){
570  iP.push_back(iP[i]);
571  zProps.push_back(zProps[i]);
572  Chi2.push_back(0);
573  trkHits.push_back(newTrkHits);
574  trkTrajs.push_back(trkTrajs[i]);
575  trkDirs.push_back(trkDirs[i]);
576  nTrkHits.push_back(nhitsnew);
577  ignoreTrk.push_back(false);
578  }
579 
580  }
581  else{
582  // if track passes check stop filtering
583  filter = false;
584  break;
585  }
586  // nhitsold = nTrkHits[i];
587  } // if last propogation
588  else{
589  // not the last propogation step stop filtering now
590  filter = false;
591  break;
592  }
593  } // end of if number of hits is the same and filter counter is greater than 0
594 
595  // update the nhits counter
596  nhitsold = nTrkHits[i];
597 
598  } // end of while loop over filtering
599  } // end if this is a track we care to filter
600  } // end of loop over tracks
601 
602  if(!initTrack){ break; }
603 
604  } // end of loop over iprop
605 
606  unsigned int loopMinHits = fMinHits;
607  if((nTimes == 0) && (fLongTrackLoop)) { //Addition by Susan Lein
608  for (unsigned int s=0; s<trkHits.size(); ++s){
609  if ((.85*nTrkHits[s]) > loopMinHits && !ignoreTrk[s]) {
610  unsigned int trkMinHits = .85*nTrkHits[s];
611  loopMinHits = trkMinHits;
612  }
613  }
614  nTimes = 1;
615  }
616  else { nTimes = 2; }
617 
618  if(trkHits.size() > 0){
619  // Find best track
620  // use the chi2/(number of hits in track) as measure of goodness of track. Track with min of this is best
621  double chiOverLength = 99999;
622  int bestTrack = -1;
623  for(unsigned int i = 0; i<trkHits.size();i++){
624  double nhits = nTrkHits[i];
625  if(nhits < loopMinHits || ignoreTrk[i]){ continue; }
626  double chiOverLengthnew = Chi2[i]/nhits;
627  if(chiOverLengthnew < chiOverLength){
628  bestTrack = i;
629  chiOverLength = chiOverLengthnew;
630  }
631  }
632  if(bestTrack >= 0){
633 
634  // make a track out of the best one
635  rb::Cluster clust(fView);
636  int firstHit = GetFirstHit(trkHits[bestTrack]);
637  rb::Track tracknew = rb::Track(clust,trkTrajs[bestTrack][firstHit][0],
638  trkTrajs[bestTrack][firstHit][1],
639  trkDirs[bestTrack][firstHit][0],
640  trkDirs[bestTrack][firstHit][1]);
641 
642 
643  // add the hits and trajectory points to the track
644  for(int ihit = 0; ihit < fNChans; ++ihit){
645  if(trkHits[bestTrack][ihit]){
646  if(zProps[bestTrack]){ tracknew.Add(fSliceChans[ihit].GetHit()); }
647  else{ tracknew.Add(fSliceChansCellSort[ihit].GetHit()); }
648  if(ihit != firstHit){ tracknew.AppendTrajectoryPoint(trkTrajs[bestTrack][ihit][0],trkTrajs[bestTrack][ihit][1]); }
649  }
650  }
651 
652  std::vector<int> extratrpts;
653  for(unsigned int itraj = 1; itraj<tracknew.NTrajectoryPoints();++itraj){
654  TVector3 oldtrpt = tracknew.TrajectoryPoint(itraj-1);
655  if(oldtrpt.X() == tracknew.TrajectoryPoint(itraj).X() &&
656  oldtrpt.Y() == tracknew.TrajectoryPoint(itraj).Y() &&
657  oldtrpt.Z() == tracknew.TrajectoryPoint(itraj).Z()){ extratrpts.push_back(itraj);}
658  }
659  std::sort(extratrpts.begin(),extratrpts.end());
660  for(int rmpt = ((int)extratrpts.size()-1);rmpt>=0;--rmpt){
661  tracknew.RemoveTrajectoryPoint(extratrpts[rmpt]);
662  }
663  if(zProps[bestTrack]){
664  // extend the track ends to the distance of closest approach of the end hits on the track.
665  // this should only be neccessary for tracks that were propagated in z, i.e. are more than one plane long
666  // get the start position and direction
667  TVector3 start = tracknew.Start();
668  TVector3 dir = tracknew.Dir();
669  double startslope = dir.X()/dir.Z();
670  if(fView == 1){startslope = dir.Y()/dir.Z(); }
671  // get the extent of the track in the start plane
672  if(startslope != 0){ // if slope is zero the fit should be correct already
673  int startCell = 10000;
674  if(startslope < 0){ startCell = -1; }
675  for(unsigned int iHit = 0; iHit<tracknew.NCell(); ++iHit){
676  if(tracknew.Cell(iHit)->Plane() == tracknew.MinPlane()){
677  if(startslope<0 && tracknew.Cell(iHit)->Cell()>startCell){startCell = tracknew.Cell(iHit)->Cell();}
678  if(startslope>0 && tracknew.Cell(iHit)->Cell()<startCell){startCell = tracknew.Cell(iHit)->Cell();}
679  }
680  }
681  if(startCell != 10000){
682  std::array<double, 2> bestStart;
683  double cpos[3];
684  const geo::CellGeo* scell = fgeom->Plane(tracknew.MinPlane())->Cell(startCell);
685  scell->GetCenter(cpos);
686  double deltaz = scell->HalfD();
687  double deltav = scell->HalfW();
688  // find the cell "best position" i.e. min dist to track or midpoint track makes through cell
689  if(fView == 0){ kgeo.BestTrackPoint(cpos[2],cpos[0],deltaz,deltav,bestStart,start.Z(),start.X(),dir.Z(),dir.X()); }
690  else if(fView == 1){ kgeo.BestTrackPoint(cpos[2],cpos[1],deltaz,deltav,bestStart,start.Z(),start.Y(),dir.Z(),dir.Y()); }
691  tracknew.SetStart(bestStart[1],bestStart[0]);
692  }
693  }
694  // get the stop position and relative direction
695  TVector3 stop = tracknew.Stop();
696  if(tracknew.NTrajectoryPoints()>1){
697  TVector3 traj = tracknew.TrajectoryPoint(tracknew.NTrajectoryPoints()-2);
698  dir = stop-traj;
699  }
700  double stopslope = dir.X()/dir.Z();
701  if(fView == 1){stopslope = dir.Y()/dir.Z(); }
702  if(stopslope != 0){ // if slope is zero the fit should be correct already
703  // get the extent of the track in the stop plane
704  int stopCell = 10000;
705  if(stopslope > 0){ stopCell = -1; }
706  for(unsigned int iHit = 0; iHit<tracknew.NCell(); ++iHit){
707  if(tracknew.Cell(iHit)->Plane() == tracknew.MaxPlane()){
708  if(stopslope<0 && tracknew.Cell(iHit)->Cell()<stopCell){stopCell = tracknew.Cell(iHit)->Cell();}
709  if(stopslope>0 && tracknew.Cell(iHit)->Cell()>stopCell){stopCell = tracknew.Cell(iHit)->Cell();}
710  }
711  }
712  if(stopCell!= 10000){
713  std::array<double, 2> bestStop;
714  double cpos[3];
715  const geo::CellGeo* stcell = fgeom->Plane(tracknew.MaxPlane())->Cell(stopCell);
716  stcell->GetCenter(cpos);
717  double deltaz = stcell->HalfD();
718  double deltav = stcell->HalfW();
719  // find the cell "best position" i.e. min dist to track or midpoint track makes through cell
720  if(fView == 0){ kgeo.BestTrackPoint(cpos[2],cpos[0],deltaz,deltav,bestStop,stop.Z(),stop.X(),dir.Z(),dir.X()); }
721  else if(fView == 1){ kgeo.BestTrackPoint(cpos[2],cpos[1],deltaz,deltav,bestStop,stop.Z(),stop.Y(),dir.Z(),dir.Y()); }
722  tracknew.RemoveTrajectoryPoint(tracknew.NTrajectoryPoints()-1);
723  tracknew.AppendTrajectoryPoint(bestStop[1],bestStop[0]);
724  }
725  }
726  }
727 
728  // if the start and stop planes are the same i.e. vertical track, pick direction so that the start is nearest to the avg pos
729  if(tracknew.MinPlane() == tracknew.MaxPlane() && tracknew.NTrajectoryPoints() > 0){
730  if(fView == 0){
731  if(abs(tracknew.Start().X() - avgX) > abs(tracknew.Stop().X() - avgX)){
732  // reverse order of trajectory points
733  std::vector<TVector3> trjpoints;
734  for(unsigned int i = 0; i<tracknew.NTrajectoryPoints(); ++i){
735  trjpoints.push_back(tracknew.TrajectoryPoint(i));
736  }
737  tracknew.ClearTrajectoryPoints();
738  reverse(trjpoints.begin(),trjpoints.end());
739  tracknew.SetStart(trjpoints[0].X(), trjpoints[0].Z());
740  for(unsigned int i = 1; i<trjpoints.size(); ++i){
741  tracknew.AppendTrajectoryPoint(trjpoints[i].X(),trjpoints[0].Z());
742  }
743  }
744  }
745  else{
746  if(abs(tracknew.Start().Y() - avgY) > abs(tracknew.Stop().Y() - avgY)){
747  // reverse order of trajectory points
748  std::vector<TVector3> trjpoints;
749  for(unsigned int i = 0; i<tracknew.NTrajectoryPoints(); ++i){
750  trjpoints.push_back(tracknew.TrajectoryPoint(i));
751  }
752  tracknew.ClearTrajectoryPoints();
753  reverse(trjpoints.begin(),trjpoints.end());
754  tracknew.SetStart(trjpoints[0].Y(), trjpoints[0].Z());
755  for(unsigned int i = 1; i<trjpoints.size(); ++i){
756  tracknew.AppendTrajectoryPoint(trjpoints[i].Y(),trjpoints[0].Z());
757  }
758  }
759  }
760  }
761 
762  // Store best track
763  // filter best track to get the correct dir vector at vertex
764  outtracks.push_back(tracknew);
765  nTimes = 0;
766 
767  // Remove track hits from signal
768  RemoveHitsFromSignal(tracknew);
769  if(fNChans < (int)fMinHits){ Done = true; }
770  fTrySingle = true;
771  }
772  else if(fTrySingle){ fTrySingle = false; }
773  else{
774  if(nTimes == 2){ Done = true; }
775  }
776  }
777  else if(fTrySingle){ fTrySingle = false; }
778  else{
779  if(nTimes == 2){ Done = true; }
780  }
781  } //end over while loop
782 
783  return outtracks;
784 
785  } // End of FindTracks function
rb::Track::NTrajectoryPoints
size_t NTrajectoryPoints() const
Definition: Track.h:83
trk::KalmanTrack::PlaneToCellSortHits
void PlaneToCellSortHits(std::vector< bool > &planeSortedHits)
Definition: KalmanTrack_module.cc:2404
rb::Track::SetStart
virtual void SetStart(TVector3 start)
Definition: Track.cxx:168
geo::CellGeo::HalfD
double HalfD() const
Definition: CellGeo.cxx:205
rb::Cluster::NCell
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
geo::CellGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
rb::Track::TrajectoryPoint
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
it
set< int >::iterator it
Definition: ProduceBadRunsByFailureMode.C:101
trk::KalmanTrack::PlaneToCellSort
void PlaneToCellSort(std::vector< std::array< double, 2 > > &planeSortedHits)
Definition: KalmanTrack_module.cc:2436
trk::KalmanTrack::FilterOnly
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)
Definition: KalmanTrack_module.cc:1755
rb::CellHit::Plane
unsigned short Plane() const
Definition: CellHit.h:39
trk::KalmanTrack::CheckTrack
void CheckTrack(KalmanGeoHelper &kgeo, std::vector< bool > &trkHits, std::vector< std::array< double, 2 > > trkTrajs, bool zProp, std::vector< bool > &newTrkHits)
Definition: KalmanTrack_module.cc:2116
geo::CellGeo::HalfW
double HalfW() const
Definition: CellGeo.cxx:191
rb::Cluster
A collection of associated CellHits.
Definition: Cluster.h:47
abs
void abs(TH1 *hist)
Definition: wrong_sign_uncertainty.C:716
rb::Track
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
geo::GeometryBase::Plane
const PlaneGeo * Plane(unsigned int i) const
Definition: GeometryBase.cxx:749
rb::Prong::Start
virtual TVector3 Start() const
Definition: Prong.h:73
site_stats_from_log.reverse
reverse
Definition: site_stats_from_log.py:192
trk::SortByCell
void SortByCell(std::vector< trk::LocatedChan > &c)
Definition: LocatedChan.cxx:78
filter
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.
Definition: TruthContain_module.cc:35
Y
Float_t Y
Definition: plot.C:38
febshutoff_auto.start
start
Definition: febshutoff_auto.py:53
rb::Cluster::Add
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
Z
Float_t Z
Definition: plot.C:38
trk::KalmanTrack::fSliceChansCellSort
std::vector< trk::LocatedChan > fSliceChansCellSort
Definition: KalmanTrack_module.cc:66
rb::CellHit::Cell
unsigned short Cell() const
Definition: CellHit.h:40
s
const XML_Char * s
Definition: expat.h:262
trk::KalmanTrack::RemoveHitsFromSignal
void RemoveHitsFromSignal(const rb::Track &track)
Definition: KalmanTrack_module.cc:1727
cet::n
std::void_t< T > n
Definition: metaprogramming.h:28
rb::Track::AppendTrajectoryPoint
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
rb::Prong::Dir
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
trk::KalmanTrack::fSliceChans
std::vector< trk::LocatedChan > fSliceChans
Definition: KalmanTrack_module.cc:65
trk::KalmanTrack::CellToPlaneSort
void CellToPlaneSort(std::vector< std::array< double, 2 > > &cellSortedHits)
Definition: KalmanTrack_module.cc:2449
rb::Track::ClearTrajectoryPoints
void ClearTrajectoryPoints()
Forget about all trajectory points.
Definition: Track.cxx:359
rb::Cluster::MinPlane
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
rb::Cluster::Cell
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
trk::KalmanTrack::RemoveSameTracks
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)
Definition: KalmanTrack_module.cc:1667
trk::KalmanTrack::IsSinglePlane
bool IsSinglePlane(const std::vector< bool > &trkHits)
Definition: KalmanTrack_module.cc:2478
dir
TDirectory * dir
Definition: macro.C:5
trk::KalmanTrack::fgeom
art::ServiceHandle< geo::Geometry > fgeom
Definition: KalmanTrack_module.cc:54
trk::SortByPlane
void SortByPlane(std::vector< trk::LocatedChan > &c)
Definition: LocatedChan.cxx:67
trk::KalmanTrack::CountHits
int CountHits(const std::vector< bool > &trkHits)
Definition: KalmanTrack_module.cc:2468
trk::KalmanTrack::GetFirstHit
int GetFirstHit(const std::vector< bool > &trkHits)
Definition: KalmanTrack_module.cc:2503
rb::Track::RemoveTrajectoryPoint
void RemoveTrajectoryPoint(unsigned int i)
Remove the ith trajectory point from the track.
Definition: Track.cxx:367
rb::Cluster::MaxPlane
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
rb::Track::Stop
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
geo::CellGeo
Encapsulate the cell geometry.
Definition: CellGeo.h:25
X
Float_t X
Definition: plot.C:38
trk::KalmanTrack::SegmentFinder
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)
Definition: KalmanTrack_module.cc:788
trk::KalmanTrack::CellToPlaneSortHits
void CellToPlaneSortHits(std::vector< bool > &cellSortedHits)
Definition: KalmanTrack_module.cc:2417
trk::KalmanTrack::FilterTracker
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)
Definition: KalmanTrack_module.cc:1104
std::array<std::vector<ProductInfo>, NumBranchTypes> const& art::ModuleBase::getConsumables ( ) const
inherited
int trk::KalmanTrack::GetFirstHit ( const std::vector< bool > &  trkHits)
private

Definition at line 2503 of file KalmanTrack_module.cc.

References fNChans, and MECModelEnuComparisons::i.

Referenced by FilterOnly(), and FindTracks().

2504  {
2505  int firstHit = -1;
2506  for(int i = 0; i < fNChans; ++i){
2507  if(trkHits[i]){
2508  firstHit = i;
2509  break;
2510  }
2511  }
2512  return firstHit;
2513  }
int trk::KalmanTrack::GetLastHit ( const std::vector< bool > &  trkHits)
private

Definition at line 2516 of file KalmanTrack_module.cc.

References DEFINE_ART_MODULE(), fNChans, and MECModelEnuComparisons::i.

2517  {
2518  int lastHit = -1;
2519  for(int i = fNChans-1; i >= 0; --i){
2520  if(trkHits[i]){
2521  lastHit = i;
2522  break;
2523  }
2524  }
2525  return lastHit;
2526  }
bool trk::KalmanTrack::IsSinglePlane ( const std::vector< bool > &  trkHits)
private

Definition at line 2478 of file KalmanTrack_module.cc.

References fNChans, fSliceChans, and MECModelEnuComparisons::i.

Referenced by FindTracks(), and SingleSegment().

2479  {
2480 
2481  int minpl = -1;
2482  int maxpl = -1;
2483  // find the first plane in this track
2484  for(int i = 0; i < fNChans; ++i){
2485  if(trkHits[i]){
2486  minpl = fSliceChans[i].Plane();
2487  break;
2488  }
2489  }
2490  // find the last plane in this track
2491  for(int i = fNChans-1; i >=0; --i){
2492  if(trkHits[i]){
2493  maxpl = fSliceChans[i].Plane();
2494  break;
2495  }
2496  }
2497  if(minpl == maxpl){ return true; }
2498  else{ return false; }
2499 
2500  }
trk::KalmanTrack::fSliceChans
std::vector< trk::LocatedChan > fSliceChans
Definition: KalmanTrack_module.cc:65
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::mayConsume ( InputTag const &  tag)
protectedinherited

Definition at line 76 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::mayConsume(), T, and getGoodRuns4SAM::tag.

77  {
78  return collector_.mayConsume<T, BT>(tag);
79  }
art::ConsumesCollector::mayConsume
ProductToken< T > mayConsume(InputTag const &)
Definition: ConsumesCollector.h:76
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:50
T
double T
Definition: Xdiff_gwt.C:5
template<typename T , BranchType BT>
void art::ModuleBase::mayConsumeMany ( )
protectedinherited

Definition at line 90 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::mayConsumeMany(), and T.

91  {
92  collector_.mayConsumeMany<T, BT>();
93  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:50
T
double T
Definition: Xdiff_gwt.C:5
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const &  tag)
inherited

Definition at line 83 of file ModuleBase.h.

References art::ModuleBase::collector_, art::ConsumesCollector::mayConsumeView(), T, and getGoodRuns4SAM::tag.

84  {
85  return collector_.mayConsumeView<T, BT>(tag);
86  }
art::ModuleBase::collector_
ConsumesCollector collector_
Definition: ModuleBase.h:50
art::ConsumesCollector::mayConsumeView
ViewToken< Element > mayConsumeView(InputTag const &)
T
double T
Definition: Xdiff_gwt.C:5
ModuleDescription const& art::ModuleBase::moduleDescription ( ) const
inherited
void trk::KalmanTrack::PlaneToCellSort ( std::vector< std::array< double, 2 > > &  planeSortedHits)
private

Definition at line 2436 of file KalmanTrack_module.cc.

References fPlaneToCell, and MECModelEnuComparisons::i.

Referenced by FindTracks().

2437  {
2438  std::vector<std::array<double, 2>> cellSortedVec(planeSortedVec.size());
2439  // loop over the planeSortedVec and set the corresponding cellSortedHit to true
2440  for(unsigned int i = 0; i < planeSortedVec.size(); ++i){
2441  // get the corresponding cell sorted hit location
2442  int cellsort = fPlaneToCell[i];
2443  cellSortedVec[cellsort] = planeSortedVec[i];
2444  }
2445  planeSortedVec = cellSortedVec;
2446  }
trk::KalmanTrack::fPlaneToCell
std::vector< int > fPlaneToCell
Definition: KalmanTrack_module.cc:67
void trk::KalmanTrack::PlaneToCellSortHits ( std::vector< bool > &  planeSortedHits)
private

Definition at line 2404 of file KalmanTrack_module.cc.

References fPlaneToCell, and MECModelEnuComparisons::i.

Referenced by FindTracks(), and SingleSegment().

2405  {
2406  std::vector<bool> cellSortedHits(planeSortedHits.size());
2407  // loop over the planeSortedHits and set the corresponding cellSortedHit to true
2408  for(unsigned int i = 0; i < planeSortedHits.size(); ++i){
2409  // get the corresponding cell sorted hit location
2410  int cellsort = fPlaneToCell[i];
2411  cellSortedHits[cellsort] = planeSortedHits[i];
2412  }
2413  planeSortedHits = cellSortedHits;
2414  }
trk::KalmanTrack::fPlaneToCell
std::vector< int > fPlaneToCell
Definition: KalmanTrack_module.cc:67
void trk::KalmanTrack::produce ( art::Event &  evt)
virtual

Reconstruction method for this module. Provides read-write access to the event data.

Implements art::EDProducer.

Definition at line 182 of file KalmanTrack_module.cc.

References avgX, avgXPos, avgXPosZRange, avgY, avgYPos, avgYPosZRange, plot_validation_datamc::c, rb::CellHit::Cell(), util::CreateAssn(), dz, fbadc, fBadChannels, fCellL, fCellW, fClusterInput, fgeom, FindTracks(), fMaxHitCut, fMaxSliceHits, fMinHits, fNChans, fObeyPreselection, fSliceChans, fView, art::DataViewImpl::getByLabel(), geo::CellGeo::GetCenter(), geo::CellGeo::HalfD(), geo::CellGeo::HalfL(), geo::CellGeo::HalfW(), MECModelEnuComparisons::i, makeTrainCVSamples::int, chaninfo::BadChanList::IsBad(), rb::IsFiltered(), it, geo::kX, geo::kY, geo::PlaneGeo::Ncells(), NCells, geo::GeometryBase::NPlanes(), rb::CellHit::Plane(), geo::GeometryBase::Plane(), art::PtrVector< T >::push_back(), art::DataViewImpl::put(), art::PtrVector< T >::reserve(), trk::LocatedChan::SetHit(), art::PtrVector< T >::size(), trk::SortByPlane(), and geo::PlaneGeo::View().

183  {
184 
185  // Declare a container for Track objects to be stored in the art::event
186  std::unique_ptr< std::vector<rb::Track> > trackcol(new std::vector<rb::Track>);
187  std::unique_ptr< art::Assns<rb::Track, rb::Cluster> > assns(new art::Assns<rb::Track, rb::Cluster>);
188 
189  // Make the helper geo function
190  KalmanGeoHelper kgeo(fBadChannels);
191 
192  // define a vector holding slices
193  art::Handle<std::vector<rb::Cluster> > slicecol;
194  evt.getByLabel(fClusterInput,slicecol);
195  //Now make it an art::PtrVector
196  art::PtrVector<rb::Cluster> slicelist;
197  for(unsigned int i = 0; i < slicecol->size(); ++i){
198  if(fObeyPreselection && rb::IsFiltered(evt,slicecol,i)){ continue; }
199  art::Ptr<rb::Cluster>slice(slicecol,i);
200  slicelist.push_back(slice);
201  }
202 
203  // get the total number of non masked off channels by looping over all planes/cells in the detector
204  int numGoodChannels = 0;
205  bool foundmaxcellx = false;
206  bool foundmaxcelly = false;
207  for(unsigned int iPlane = 0; iPlane<fgeom->NPlanes();++iPlane){
208  const geo::PlaneGeo* p = fgeom->Plane(iPlane);
209  if(!foundmaxcellx && p->View() == geo::kX){
210  NCells[0] = p->Ncells()-1;
211  foundmaxcellx = true;
212  }
213  else if(!foundmaxcelly && p->View() == geo::kY){
214  NCells[1] = p->Ncells()-1;
215  foundmaxcelly = true;
216  }
217  for(unsigned int iCell = 0;iCell<p->Ncells();++iCell){
218  if(!(fbadc->IsBad(iPlane,iCell))){ ++numGoodChannels; }
219  }
220  }
221 
222  double xCellW = 0.0;
223  double xCellL = 0.0;
224  double yCellW = 0.0;
225  double yCellL = 0.0;
226 
227  // Loop over slices to try to fit tracks in each individual slice
228  // ignoring the noise slice
229  for(unsigned int iSlice = 0; iSlice < slicelist.size(); iSlice++){
230  // skip noise slices and big slices
231  if(slicelist[iSlice]->IsNoise() || slicelist[iSlice]->NCell() > fMaxHitCut){ continue; }
232 
233  // Check if this slice is a likely air shower. If so, do not attempt reco
234  if(slicelist[iSlice]->NCell() < numGoodChannels*fMaxSliceHits){
235  avgXPosZRange.resize(0);
236  avgXPos.resize(0);
237  avgYPosZRange.resize(0);
238  avgYPos.resize(0);
239  avgX = slicelist[iSlice]->MeanX();
240  avgY = slicelist[iSlice]->MeanY();
241 
242  int minplanex = 10000;
243  int maxplanex = 0;
244 
245  // Get the x cell hits out of the slice
246  std::vector<trk::LocatedChan> xcellchans;
247  xcellchans.reserve(slicelist[iSlice]->NXCell());
248  double zxmin = 100000;
249  double zxmax = 0;
250  for(unsigned int nx = 0; nx < slicelist[iSlice]->NXCell(); ++nx){
251  art::Ptr<rb::CellHit> xhit = slicelist[iSlice]->XCell(nx);
252  LocatedChan xchan(xhit->Plane(),xhit->Cell());
253  xchan.SetHit(xhit);
254  double xyz[3];
255  const geo::CellGeo* c = fgeom->Plane(xchan.Plane())->Cell(xchan.Cell());
256  c->GetCenter(xyz);
257  xchan.SetCenter(xyz[0],xyz[1],xyz[2]);
258  xchan.SetHalfWidths(c->HalfW(),c->HalfL(),c->HalfD());
259  xcellchans.push_back(xchan);
260  if(xchan.Plane() < minplanex){ minplanex = xchan.Plane(); }
261  if(xchan.Plane() > maxplanex){ maxplanex = xchan.Plane(); }
262  if(nx == 0){
263  xCellW = c->HalfW();
264  xCellL = c->HalfD();
265  }
266  if(xyz[2]>zxmax){ zxmax = xyz[2]; }
267  if(xyz[2]<zxmin){ zxmin = xyz[2]; }
268  }
269 
270  // get a vector that holds the average x position for every 100 cm in z
271  trk::SortByPlane(xcellchans);
272  // how many break points in z
273  unsigned int nxzbreak = zxmax-zxmin > 0? (int)floor((zxmax-zxmin)/100) : 0;
274  // decide where to break the z positions up
275  double dz = (zxmax-zxmin)/((double)nxzbreak);
276  double xzbreaks[nxzbreak+2];
277  xzbreaks[0] = zxmin;
278  for(unsigned int ibrk = 1; ibrk<=nxzbreak; ++ibrk){
279  xzbreaks[ibrk] = zxmin+((double)ibrk)*dz;
280  }
281  xzbreaks[nxzbreak+1] = zxmax;
282  for(unsigned int ibrk = 1; ibrk < nxzbreak+2; ++ibrk){
283  double xtot = 0;
284  double ntot = 0;
285  for(unsigned int xchan = 0;xchan<xcellchans.size();++xchan){
286  double xyz[3];
287  xcellchans[xchan].GetCenter(xyz);
288  if(xyz[2] >= xzbreaks[ibrk-1] && xyz[2] <=xzbreaks[ibrk]){
289  xtot+=xyz[0];
290  ntot+=1.0;
291  }
292  }
293  std::vector<double> zrng;
294  zrng.push_back(xzbreaks[ibrk-1]);
295  zrng.push_back(xzbreaks[ibrk]);
296  avgXPosZRange.push_back(zrng);
297  if(ntot>0){ avgXPos.push_back(xtot/ntot); }
298  else{ avgXPos.push_back(avgX); }
299  }
300 
301  int minplaney = 10000;
302  int maxplaney = 0;
303 
304  // Get the y cell hits out of the slice
305  std::vector<trk::LocatedChan> ycellchans;
306  double zymin = 100000;
307  double zymax = 0;
308  for(unsigned int ny = 0; ny < slicelist[iSlice]->NYCell(); ++ny){
309  art::Ptr<rb::CellHit> yhit = slicelist[iSlice]->YCell(ny);
310  LocatedChan ychan(yhit->Plane(),yhit->Cell());
311  ychan.SetHit(yhit);
312  double xyz[3];
313  const geo::CellGeo* c = fgeom->Plane(ychan.Plane())->Cell(ychan.Cell());
314  c->GetCenter(xyz);
315  ychan.SetCenter(xyz[0],xyz[1],xyz[2]);
316  ychan.SetHalfWidths(c->HalfL(),c->HalfW(),c->HalfD());
317  ycellchans.push_back(ychan);
318  if(ychan.Plane() < minplaney){ minplaney = ychan.Plane(); }
319  if(ychan.Plane() > maxplaney){ maxplaney = ychan.Plane(); }
320  if(ny == 0){
321  yCellW = c->HalfW();
322  yCellL = c->HalfD();
323  }
324  if(xyz[2]>zymax){ zymax = xyz[2]; }
325  if(xyz[2]<zymin){ zymin = xyz[2]; }
326  }
327 
328  // how many break points in yz
329  int nyzbreak = zymax-zymin > 0? (int)floor((zymax-zymin)/100) : 0;
330  // decide where to break the z positions up
331  dz = (zymax-zymin)/((double)nyzbreak);
332  std::vector<double> yzbreaks;
333  yzbreaks.push_back(zymin);
334  for(int ibrk = 1; ibrk<=nyzbreak; ++ibrk){
335  yzbreaks.push_back(zymin+((double)ibrk)*dz);
336  }
337  yzbreaks.push_back(zymax);
338  for(unsigned int ibrk = 1; ibrk < yzbreaks.size(); ++ibrk){
339  double ytot = 0;
340  double ntot = 0;
341  for(unsigned int ychan = 0;ychan<ycellchans.size();++ychan){
342  double xyz[3];
343  ycellchans[ychan].GetCenter(xyz);
344  if(xyz[2] >= yzbreaks[ibrk-1] && xyz[2] <=yzbreaks[ibrk]){
345  ytot+=xyz[1];
346  ntot+=1.0;
347  }
348  }
349  std::vector<double> zrng;
350  zrng.push_back(yzbreaks[ibrk-1]);
351  zrng.push_back(yzbreaks[ibrk]);
352  avgYPosZRange.push_back(zrng);
353  if(ntot>0){ avgYPos.push_back(ytot/ntot); }
354  else{ avgYPos.push_back(avgY); }
355  }
356 
357  // Define the vectors to hold the track info from the new FindTracks function
358  std::vector<rb::Track> xTracks;
359  fSliceChans = xcellchans;
360  fNChans = fSliceChans.size();
361  fView = geo::kX;
362  fCellW = xCellW;
363  fCellL = xCellL;
364 
365  if(xcellchans.size() >= fMinHits){ xTracks = FindTracks(kgeo, xcellchans); }
366 
367  // Define the vectors to hold the track info from the new FindTracks function
368  std::vector<rb::Track> yTracks;
369  fSliceChans = ycellchans;
370  fNChans = fSliceChans.size();
371  fView = geo::kY;
372  fCellW = yCellW;
373  fCellL = yCellL;
374 
375  if(ycellchans.size() >= fMinHits){yTracks = FindTracks(kgeo, ycellchans); }
376 
377  std::vector<rb::Track> alltracks;
378  for(unsigned int i = 0; i<xTracks.size();++i){
379  xTracks[i].SetID(iSlice);
380  alltracks.push_back(xTracks[i]);
381  }
382  for(unsigned int i = 0; i<yTracks.size();++i){
383  yTracks[i].SetID(iSlice);
384  alltracks.push_back(yTracks[i]);
385  }
386  for(std::vector<rb::Track>::const_iterator it = alltracks.begin(); it != alltracks.end(); ++it){
387  trackcol->push_back(*it);
388  util::CreateAssn(evt,*trackcol,slicelist[iSlice],*assns);
389  }
390  } // end of if statement on requirement on the number of hits in slice
391  } // end of loop over slices
392 
393  evt.put(std::move(trackcol));
394  evt.put(std::move(assns));
395 
396  }
trk::KalmanTrack::avgXPosZRange
std::vector< std::vector< double > > avgXPosZRange
Definition: KalmanTrack_module.cc:57
art::PtrVector::reserve
void reserve(size_type n)
Definition: PtrVector.h:337
geo::CellGeo::HalfL
double HalfL() const
Definition: CellGeo.cxx:198
util::CreateAssn
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.
geo::CellGeo::HalfD
double HalfD() const
Definition: CellGeo.cxx:205
geo::CellGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
it
set< int >::iterator it
Definition: ProduceBadRunsByFailureMode.C:101
rb::CellHit::Plane
unsigned short Plane() const
Definition: CellHit.h:39
trk::KalmanTrack::FindTracks
std::vector< rb::Track > FindTracks(KalmanGeoHelper &kgeo, std::vector< trk::LocatedChan > sliceChans)
Definition: KalmanTrack_module.cc:399
p
const char * p
Definition: xmltok.h:285
trk::KalmanTrack::fClusterInput
std::string fClusterInput
Definition: KalmanTrack_module.cc:44
geo::kX
Vertical planes which measure X.
Definition: PlaneGeo.h:28
geo::PlaneGeo::Ncells
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
geo::CellGeo::HalfW
double HalfW() const
Definition: CellGeo.cxx:191
trk::KalmanTrack::avgXPos
std::vector< double > avgXPos
Definition: KalmanTrack_module.cc:56
geo::GeometryBase::Plane
const PlaneGeo * Plane(unsigned int i) const
Definition: GeometryBase.cxx:749
geo::kY
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
rb::CellHit::Cell
unsigned short Cell() const
Definition: CellHit.h:40
geo::PlaneGeo::View
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
dz
double dz[NP][NC]
Definition: checkMeasured_Target.C:20
art::PtrVector::push_back
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
trk::KalmanTrack::fbadc
art::ServiceHandle< chaninfo::BadChanList > fbadc
Definition: KalmanTrack_module.cc:55
geo::PlaneGeo
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
plot_validation_datamc.c
c
Definition: plot_validation_datamc.py:198
trk::KalmanTrack::fSliceChans
std::vector< trk::LocatedChan > fSliceChans
Definition: KalmanTrack_module.cc:65
rb::IsFiltered
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
art::PtrVector::size
size_type size() const
Definition: PtrVector.h:302
trk::KalmanTrack::avgYPos
std::vector< double > avgYPos
Definition: KalmanTrack_module.cc:58
trk::KalmanTrack::avgYPosZRange
std::vector< std::vector< double > > avgYPosZRange
Definition: KalmanTrack_module.cc:59
trk::KalmanTrack::fgeom
art::ServiceHandle< geo::Geometry > fgeom
Definition: KalmanTrack_module.cc:54
trk::SortByPlane
void SortByPlane(std::vector< trk::LocatedChan > &c)
Definition: LocatedChan.cxx:67
geo::GeometryBase::NPlanes
unsigned int NPlanes() const
Definition: GeometryBase.cxx:460
geo::CellGeo
Encapsulate the cell geometry.
Definition: CellGeo.h:25
art::DataViewImpl::put
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
chaninfo::BadChanList::IsBad
bool IsBad(int plane, int cell)
Definition: BadChanList_service.cc:563
void trk::KalmanTrack::reconfigure ( const fhicl::ParameterSet &  p)

Definition at line 162 of file KalmanTrack_module.cc.

References fBadChannels, fClusterInput, fDeltaChiAllowed, fGapAllowed, fLargeSliceHits, fLongTrackLoop, fMaxHitCut, fMaxSliceHits, fMinGapProb, fMinHits, fObeyPreselection, fWhoaNelly, fhicl::ParameterSet::get(), and string.

Referenced by KalmanTrack().

163  {
164  fMaxSliceHits = pset.get< double >("MaxSliceHits") ;
165  fMaxHitCut = pset.get< double >("MaxHitCut") ;
166  fWhoaNelly = pset.get< int >("WhoaNelly") ;
167  fClusterInput = pset.get< std::string >("ClusterInput") ;
168  fDeltaChiAllowed = pset.get< double >("DeltaChiAllowed") ;
169  fMinHits = pset.get< int >("MinHits") ;
170  fGapAllowed = pset.get< int >("GapAllowed") ;
171  fLargeSliceHits = pset.get< int >("LargeSliceHits") ;
172  fBadChannels = pset.get< bool >("BadChannels") ;
173  fLongTrackLoop = pset.get< bool >("LongTrackLoop") ;
174  fMinGapProb = pset.get< double >("MinGapProb") ;
175  fObeyPreselection = pset.get< bool >("ObeyPreselection");
176  }
trk::KalmanTrack::fClusterInput
std::string fClusterInput
Definition: KalmanTrack_module.cc:44
string
enum BeamMode string
void trk::KalmanTrack::RemoveHitsFromSignal ( const rb::Track &  track)
private

Following is added for new hit storage method

Definition at line 1727 of file KalmanTrack_module.cc.

References rb::Cluster::Cell(), fNChans, fSliceChans, MECModelEnuComparisons::i, calib::j, and rb::Cluster::NCell().

Referenced by FindTracks().

1728  {
1729 
1730  /// Following is added for new hit storage method
1731  std::vector<trk::LocatedChan> newfSignalChan;
1732  // loop over all the signal hits
1733  for(int i = 0; i<fNChans; ++i){
1734  // get current hit
1735  art::Ptr<rb::CellHit> chanhit = fSliceChans[i].GetHit();
1736  // define a boolean to track if signal hit is in track hit
1737  bool inTrack = false;
1738  for(unsigned int j = 0; j < track.NCell(); ++j){
1739  // get track hit
1740  art::Ptr<rb::CellHit> thit = track.Cell(j);
1741  if(chanhit == thit){
1742  inTrack = true;
1743  break;
1744  }
1745  }
1746  if(!inTrack){ newfSignalChan.push_back(fSliceChans[i]); }
1747  }
1748  fNChans = newfSignalChan.size();
1749  fSliceChans.resize(fNChans);
1750  fSliceChans = newfSignalChan;
1751 
1752  }
rb::Cluster::NCell
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
trk::KalmanTrack::fSliceChans
std::vector< trk::LocatedChan > fSliceChans
Definition: KalmanTrack_module.cc:65
calib::j
const double j
Definition: BetheBloch.cxx:29
rb::Cluster::Cell
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
void trk::KalmanTrack::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 
)
private

Definition at line 1667 of file KalmanTrack_module.cc.

References fNChans.

Referenced by FindTracks().

1672  {
1673  int nTracks = trkHits.size();
1674  if(nTracks == 0) return;
1675  // set up a vector to track whether or not an input track is the same as another track.
1676  bool sameTrack[nTracks];
1677  for(int st = 0; st < nTracks; ++st) sameTrack[st] = false;
1678 
1679  int nSame = 0;
1680  for(int iTrack = 0; iTrack < nTracks; iTrack++){
1681  //define a variable to track whether this track is the same as another or not.
1682  bool inTrackDone = false;
1683  if(sameTrack[iTrack] || ignoreTrk[iTrack]){inTrackDone = true; }
1684  // set a counter of the track to compare the input track to.
1685  int cTrack = iTrack+1;
1686  // loop over all other tracks
1687  while(cTrack < nTracks && !inTrackDone){
1688  if(!sameTrack[cTrack] && !ignoreTrk[cTrack]){
1689  // assume these are the same tracks until proven otherwise
1690  bool same = true;
1691  // first check if these are both vertical or horizontally propogated
1692  if(zProps[iTrack] == zProps[cTrack] && nTrkHits[iTrack] == nTrkHits[cTrack]){
1693  // compare if the same trkHits
1694  for(int ihit = 0; ihit < fNChans; ++ihit){
1695  if(trkHits[iTrack][ihit] != trkHits[cTrack][ihit]){
1696  same = false;
1697  break;
1698  }
1699  }
1700  }
1701  else{ same = false; }
1702  // if the same track keep the one with the smaller chi2
1703  if(same){
1704  if(chi2[iTrack] > chi2[cTrack]){
1705  sameTrack[iTrack] = true;
1706  ignoreTrk[iTrack] = true;
1707  // move onto next input track if this track is to be deleted
1708  inTrackDone = true;
1709  }
1710  else{
1711  sameTrack[cTrack] = true;
1712  ignoreTrk[cTrack] = true;
1713  }
1714  }
1715 
1716  }
1717  cTrack++; // move on to next track for comparison
1718  //end loop over the comparison tracks
1719  }
1720  if(sameTrack[iTrack]){ ++nSame; }
1721  //end loop over track to remove
1722  }
1723 
1724  }
chi2
double chi2()
Definition: cana_time_peak.C:532
int trk::KalmanTrack::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 
)
private

Definition at line 788 of file KalmanTrack_module.cc.

References trk::KalmanGeoHelper::CountMissedCellsOnLine(), delta, makeTrainCVSamples::dirs, fCellL, fCellW, fGapAllowed, fNChans, fPlaneToCell, fSliceChans, fTrySingle, fView, fWhoaNelly, hits(), m, P, SingleSegment(), and std::sqrt().

Referenced by FindTracks().

794  {
795  int ntracks = 0;
796  int maxgap = fGapAllowed;
797  std::vector<bool> hitsTemplate(fNChans,false);
798  std::vector<std::array<double, 2> > trajs(fNChans);
799  std::vector<std::array<double, 2> > dirs(fNChans);
800  std::vector<bool> usedVert(fNChans,false);
801  std::vector<bool> usedHor(fNChans,false);
802 
803  int firstiHitPlane = fSliceChans[fNChans-1].Plane();
804 
805  if(fTrySingle){
806  int firsthitplane = 0;
807  SingleSegment(iP,zProp,trkHits,trkTrajs,trkDirs,firsthitplane);
808  if(iP.size() > 0){
809  firstiHitPlane = 0;
810  ++ntracks;
811  }
812 
813  int planeextent = fSliceChans[fNChans-1].Plane() - fSliceChans[0].Plane();
814 
815  // Check if also grouping just the last half of the cells works better
816  firsthitplane = planeextent/2+fSliceChans[0].Plane();
817  SingleSegment(iP,zProp,trkHits,trkTrajs,trkDirs,firsthitplane);
818  if((int)iP.size() > ntracks){
819  ++ntracks;
820  // if we didn't find a usable track using all the hits update firstiHitPlane here so we can consider
821  if(ntracks < 2){ firstiHitPlane = firsthitplane; }
822 
823  }
824  if(ntracks < 2){
825  // try another division of hits
826  firsthitplane = firsthitplane + planeextent/4;
827  SingleSegment(iP,zProp,trkHits,trkTrajs,trkDirs,firsthitplane);
828  if((int)iP.size() > ntracks){
829  ++ntracks;
830  // if we didn't find a usable track using all the hits update firstiHitPlane here so we can consider
831  if(ntracks < 2){ firstiHitPlane = firsthitplane; }
832  }
833  }
834 
835 
836  if(ntracks > 1){ firstiHitPlane = 0; }
837 
838  }
839 
840  // If we didn't find any tracks this way don't try it anymore
841  if(ntracks == 0){ fTrySingle = false; }
842 
843  // loop over all points cell hits and consider them as the starting point
844  for(int iHit = fNChans-1; iHit > 0; --iHit){
845  // get this hits location
846  int iPlane = fSliceChans[iHit].Plane();
847  if(fSliceChans[iHit].Plane() >= firstiHitPlane){ continue; }
848  int iCell = fSliceChans[iHit].Cell();
849  // loop over all other cell hits and consider them as a compatible point
850  for(int fHit = iHit-1; fHit >= 0; --fHit){
851 
852  // check that number of seed is not less than maximum;
853  if(ntracks >= fWhoaNelly){ return ntracks; }
854 
855  //make loose cut on number of planes between the hits
856  int fPlane = fSliceChans[fHit].Plane();
857  int deltap = (iPlane-fPlane-2)/2;
858 
859  // if too large a gap in planes skip
860  if(deltap > maxgap){ continue; }
861 
862  // get information about the cell hits
863  int fCell = fSliceChans[fHit].Cell();
864  //make vectors to hold position information.
865  double ixyz[3]; // xyz[0] = x, xyz[1] = y, xyz[2] = z
866  fSliceChans[iHit].GetCenter(ixyz);
867  double fxyz[3];
868  fSliceChans[fHit].GetCenter(fxyz);
869 
870  // check if we can make a track out of these two hits
871  // Break up into 3 cases, completely veritcal in view, completely horizontal in z, and else
872  // First two cases can take some short cuts
873  if(iPlane == fPlane && !usedVert[iHit]){
874  // hits on the same plane check number of cells between hits
875  if(TMath::Abs(iCell-fCell)<=1){
876  // make a track segment out of these hits
877  // hits on same plane want to propogate by cell/not plane
878  zProp.push_back(false);
879  double delta = TMath::Abs(ixyz[fView] - fxyz[fView]);
880  std::vector<bool> hits = hitsTemplate;
881  int celliHitid = fPlaneToCell[iHit];
882  int cellfHitid = fPlaneToCell[fHit];
883  hits[celliHitid] = true;
884  hits[cellfHitid] = true;
885  trkHits.push_back(hits);
886  trajs[celliHitid][0] = ixyz[fView];
887  trajs[celliHitid][1] = ixyz[2];
888  trajs[cellfHitid][0] = fxyz[fView];
889  trajs[cellfHitid][1] = fxyz[2];
890  trkTrajs.push_back(trajs);
891  dirs[celliHitid][0] = 1.0;
892  dirs[celliHitid][1] = 0.0;
893  dirs[cellfHitid] = dirs[celliHitid];
894  trkDirs.push_back(dirs);
895  double Pconst = fCellL*fCellL/3.0;
896  std::vector<double> P(3);
897  P[0] = Pconst;
898  P[1] = Pconst/delta;
899  P[2] = 2*P[1]/delta;
900  iP.push_back(P);
901  usedVert[iHit] = true;
902  ++ntracks;
903  }
904  }
905  else if(iCell == fCell && !usedHor[iHit]){
906  // hits in same cell number different plane
907  // Already checked distance between plane, no more checks to make a track segment
908  // hits exactly horizontal to each other in view so propogate by plane
909  zProp.push_back(true);
910  double deltaZ = (ixyz[2]-fCellL)-(fxyz[2]+fCellL);
911  std::vector<double> P(3);
912  double Pconst = fCellW*fCellW/3.0;
913  P[0] = Pconst;
914  P[1] = Pconst/deltaZ;
915  P[2] = 2*P[1]/deltaZ;
916  iP.push_back(P);
917  usedHor[iHit] = true;
918  std::vector<bool> hits = hitsTemplate;
919  hits[iHit] = true;
920  hits[fHit] = true;
921  trkHits.push_back(hits);
922  trajs[iHit][0] = ixyz[fView];
923  trajs[iHit][1] = ixyz[2];
924  trajs[fHit][0] = fxyz[fView];
925  trajs[fHit][1] = fxyz[2];
926  trkTrajs.push_back(trajs);
927  dirs[iHit][0] = 0.0;
928  dirs[iHit][1] = 1.0;
929  dirs[fHit] = dirs[iHit];
930  trkDirs.push_back(dirs);
931  ++ntracks;
932  }
933  else if(iPlane != fPlane && iCell != fCell){
934 
935  double xyz[3] = {fxyz[0], fxyz[1], fxyz[2]+fCellL-0.001};
936  double xyz1[3] = {ixyz[0], ixyz[1], ixyz[2]-fCellL-0.001};
937  double deltaZ = ixyz[2]-fxyz[2]-2.0*fCellL;
938  double delta = ixyz[fView]-fxyz[fView];
939  double m = (1.0/(deltaZ+4*fCellL));
940  // from these potentially matched hits calculate how many hits lie between them
941  if(iCell>fCell){
942  xyz[fView]+=(fCellW+0.001);
943  xyz1[fView]-=(fCellW+0.001);
944  delta-=(2.0*fCellW);
945  m*=TMath::Abs((delta-4*fCellW));
946  }
947  else if(iCell<fCell){
948  xyz[fView]-=(fCellW+0.001);
949  xyz1[fView]+=(fCellW+0.001);
950  delta+=(2.0*fCellW);
951  m*=TMath::Abs((delta+4*fCellW));
952  }
953 
954  // Check how many cells the track passes through between iHit and fHit
955  // By construction adjacent plane hits miss zero times
956  int numMissed = 0;
957  if(iPlane>(2+fPlane)){
958  numMissed = kgeo.CountMissedCellsOnLine(fView,xyz[fView],xyz[2],xyz1[fView],xyz1[2],maxgap);
959  }
960  if(numMissed <= maxgap){
961  double slope = delta/deltaZ;
962  double dxyz[3] = {0, 0, 1/sqrt(1+slope*slope)};
963  dxyz[fView] = slope/sqrt(1+slope*slope);
964 
965  std::vector<double> P(3);
966  double Pconst = (fCellL*m+fCellW)*(fCellL*m+fCellW)/3;
967  P[0] = Pconst;
968  P[1] = Pconst/deltaZ;
969  P[2] = 2*P[1]/deltaZ;
970  std::vector<bool> hits = hitsTemplate;
971  hits[iHit] = true;
972  hits[fHit] = true;
973  trajs[iHit][0] = ixyz[fView];
974  trajs[iHit][1] = ixyz[2];
975  trajs[fHit][0] = fxyz[fView];
976  trajs[fHit][1] = fxyz[2];
977  dirs[iHit][0] = dxyz[fView];
978  dirs[iHit][1] = dxyz[2];
979  dirs[fHit] = dirs[iHit];
980 
981  trkHits.push_back(hits);
982  zProp.push_back(true);
983  iP.push_back(P);
984  trkTrajs.push_back(trajs);
985  trkDirs.push_back(dirs);
986  ++ntracks;
987  }
988  }
989  } // end of loop over fHit
990  } // end of loop over iHit
991  return ntracks;
992 
993  } // end of SegmentFinder
delta
double delta
Definition: runWimpSim.h:98
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
m
#define m
Definition: MakeGENIEReweightTable_module.cc:421
P
#define P(a, b, c, d, e, x)
hits
void hits()
Definition: readHits.C:15
trk::KalmanTrack::fSliceChans
std::vector< trk::LocatedChan > fSliceChans
Definition: KalmanTrack_module.cc:65
trk::KalmanTrack::SingleSegment
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)
Definition: KalmanTrack_module.cc:997
trk::KalmanTrack::fPlaneToCell
std::vector< int > fPlaneToCell
Definition: KalmanTrack_module.cc:67
void art::ModuleBase::setModuleDescription ( ModuleDescription const &  )
inherited
double trk::KalmanTrack::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 
)
private

Definition at line 997 of file KalmanTrack_module.cc.

References channels, makeTrainCVSamples::dirs, fNChans, fSliceChans, fSliceChansCellSort, fView, hits(), MECModelEnuComparisons::i, IsSinglePlane(), geo::LinFitMinDperp(), P, PlaneToCellSortHits(), w, and zz.

Referenced by SegmentFinder().

1003  {
1004 
1005  std::vector<trk::LocatedChan> channels = fSliceChans;
1006 
1007  // loop over all the hits and mark the indicated ones as hits to include on the track
1008  // also grab the view and z positions as well as the pe to get an initial fit
1009  std::vector<double> vv;
1010  std::vector<double> zz;
1011  std::vector<double> w;
1012  double petot = 0.0;
1013  std::vector<bool> hits(fNChans,false);
1014  // get the max/min view and z while we're at it
1015  double minz = 100000;
1016  double maxz = 0;
1017  double minv = 10000;
1018  double maxv = -10000;
1019  for(int ihit = 0; ihit < fNChans; ++ihit){
1020  if(channels[ihit].Plane() >= firsthitplane){
1021  hits[ihit] = true;
1022  double hitxyz[3];
1023  channels[ihit].GetCenter(hitxyz);
1024  vv.push_back(hitxyz[fView]);
1025  zz.push_back(hitxyz[2]);
1026  w.push_back(channels[ihit].GetHit()->PE());
1027  petot+=channels[ihit].GetHit()->PE();
1028  if(hitxyz[fView] > maxv){ maxv = hitxyz[fView]; }
1029  else if (hitxyz[fView] < minv){ minv = hitxyz[fView]; }
1030  if(hitxyz[2] > maxz){ maxz = hitxyz[2]; }
1031  else if(hitxyz[2] < minz){ minz = hitxyz[2]; }
1032  }
1033  }
1034 
1035  // do the fit
1036  double linfitchi2 = 0.0;
1037  int nhits = vv.size();
1038  if(nhits < 3){ return 0.0; }
1039  for(int i = 0; i < nhits; ++i){
1040  w[i] = w[i]/petot;
1041  }
1042  double z1[1],v1[1],z2[1],v2[1];
1043  linfitchi2 =geo::LinFitMinDperp(zz,vv,w,z1,v1,z2,v2);
1044 
1045  double chindf = linfitchi2/(nhits-2);
1046  // is this fit good enough to consider this a track segement?
1047  if(chindf > 3){ return chindf; }
1048 
1049  // Are all the hits on a single plane?
1050  bool prop = IsSinglePlane(hits);
1051  if(!prop){
1052  // sort the hits by cells
1053  channels = fSliceChansCellSort;
1054  PlaneToCellSortHits(hits);
1055  }
1056 
1057  double deltav = v2[0] - v1[0];
1058  double deltaz = z2[0] - z1[0];
1059  double dirv = deltav/(deltav*deltav+deltaz*deltaz);
1060  double dirz = deltaz/(deltav*deltav+deltaz*deltaz);
1061 
1062  // set the trajectory points and directions of the track
1063  std::vector<std::array<double, 2> > trajs(fNChans);
1064  std::vector<std::array<double, 2> > dirs(fNChans);
1065 
1066  for(int ihit = 0; ihit < fNChans; ++ihit){
1067 
1068  if(hits[ihit]){
1069  double hitxyz[3];
1070  channels[ihit].GetCenter(hitxyz);
1071  if(prop){
1072  trajs[ihit][1] = hitxyz[2];
1073  trajs[ihit][0] = (dirv/dirz)*(hitxyz[2]-z1[0])+v1[0];
1074  }
1075  else{
1076  trajs[ihit][0] = hitxyz[fView];
1077  trajs[ihit][1] = (dirz/dirv)*(hitxyz[fView]-v1[0])+z1[0];
1078  }
1079  dirs[ihit][0] = dirv;
1080  dirs[ihit][1] = dirz;
1081  }
1082 
1083  }
1084 
1085  double vdelta = maxv-minv;
1086  double zdelta = maxz-minz;
1087  // set the initial covariance matrix
1088  double Pconst = vdelta*vdelta/12.0;
1089  std::vector<double> P(3);
1090  P[0] = Pconst;
1091  P[1] = Pconst/zdelta;
1092  P[2] = 2*P[1]/zdelta;
1093 
1094  trkHits.push_back(hits);
1095  trkTrajs.push_back(trajs);
1096  trkDirs.push_back(dirs);
1097  zProp.push_back(prop);
1098  iP.push_back(P);
1099  return chindf;
1100  } // end of SingleSegments
trk::KalmanTrack::PlaneToCellSortHits
void PlaneToCellSortHits(std::vector< bool > &planeSortedHits)
Definition: KalmanTrack_module.cc:2404
zz
Double_t zz
Definition: plot.C:277
trk::KalmanTrack::fSliceChansCellSort
std::vector< trk::LocatedChan > fSliceChansCellSort
Definition: KalmanTrack_module.cc:66
P
#define P(a, b, c, d, e, x)
channels
std::map< ToFCounter, std::vector< unsigned int > > channels
Definition: ToFPositionRecoAnalysis_module.cc:81
hits
void hits()
Definition: readHits.C:15
geo::LinFitMinDperp
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...
Definition: Geo.cxx:449
trk::KalmanTrack::fSliceChans
std::vector< trk::LocatedChan > fSliceChans
Definition: KalmanTrack_module.cc:65
trk::KalmanTrack::IsSinglePlane
bool IsSinglePlane(const std::vector< bool > &trkHits)
Definition: KalmanTrack_module.cc:2478
w
Float_t w
Definition: plot.C:20
void art::ModuleBase::sortConsumables ( std::string const &  current_process_name)
inherited
std::string art::EDProducer::workerType ( ) const
inherited

Referenced by art::EDProducer::EDProducer().

Member Data Documentation

double trk::KalmanTrack::avgX
private

Definition at line 60 of file KalmanTrack_module.cc.

Referenced by FindTracks(), and produce().

std::vector<double> trk::KalmanTrack::avgXPos
private

Definition at line 56 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), and produce().

std::vector<std::vector<double> > trk::KalmanTrack::avgXPosZRange
private

Definition at line 57 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), and produce().

double trk::KalmanTrack::avgY
private

Definition at line 61 of file KalmanTrack_module.cc.

Referenced by FindTracks(), and produce().

std::vector<double> trk::KalmanTrack::avgYPos
private

Definition at line 58 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), and produce().

std::vector<std::vector<double> > trk::KalmanTrack::avgYPosZRange
private

Definition at line 59 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), and produce().

art::ServiceHandle<chaninfo::BadChanList> trk::KalmanTrack::fbadc
private

Definition at line 55 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), FilterTracker(), and produce().

bool trk::KalmanTrack::fBadChannels
private

Definition at line 49 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), FilterTracker(), produce(), and reconfigure().

double trk::KalmanTrack::fCellL
private

Definition at line 70 of file KalmanTrack_module.cc.

Referenced by FilterOnly(), FilterTracker(), produce(), and SegmentFinder().

double trk::KalmanTrack::fCellW
private

Definition at line 69 of file KalmanTrack_module.cc.

Referenced by FilterOnly(), FilterTracker(), produce(), and SegmentFinder().

std::string trk::KalmanTrack::fClusterInput
private

Definition at line 44 of file KalmanTrack_module.cc.

Referenced by produce(), and reconfigure().

double trk::KalmanTrack::fDeltaChiAllowed
private

Definition at line 45 of file KalmanTrack_module.cc.

Referenced by FilterOnly(), FilterTracker(), and reconfigure().

int trk::KalmanTrack::fGapAllowed
private

Definition at line 47 of file KalmanTrack_module.cc.

Referenced by FilterTracker(), reconfigure(), and SegmentFinder().

art::ServiceHandle<geo::Geometry> trk::KalmanTrack::fgeom
private

Definition at line 54 of file KalmanTrack_module.cc.

Referenced by FilterTracker(), FindTracks(), and produce().

int trk::KalmanTrack::fLargeSliceHits
private

Definition at line 48 of file KalmanTrack_module.cc.

Referenced by FindTracks(), and reconfigure().

bool trk::KalmanTrack::fLongTrackLoop
private

Definition at line 50 of file KalmanTrack_module.cc.

Referenced by FindTracks(), and reconfigure().

double trk::KalmanTrack::fMaxHitCut
private

Definition at line 42 of file KalmanTrack_module.cc.

Referenced by produce(), and reconfigure().

double trk::KalmanTrack::fMaxSliceHits
private

Definition at line 41 of file KalmanTrack_module.cc.

Referenced by produce(), and reconfigure().

double trk::KalmanTrack::fMinGapProb
private

Definition at line 51 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), and reconfigure().

unsigned int trk::KalmanTrack::fMinHits
private

Definition at line 46 of file KalmanTrack_module.cc.

Referenced by FindTracks(), produce(), and reconfigure().

int trk::KalmanTrack::fNChans
private

Definition at line 64 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), CountHits(), CreateHitMaps(), FilterOnly(), FilterTracker(), FindTracks(), GetFirstHit(), GetLastHit(), IsSinglePlane(), produce(), RemoveHitsFromSignal(), RemoveSameTracks(), SegmentFinder(), and SingleSegment().

bool trk::KalmanTrack::fObeyPreselection
private

Definition at line 52 of file KalmanTrack_module.cc.

Referenced by produce(), and reconfigure().

std::vector<int> trk::KalmanTrack::fPlaneToCell
private

Definition at line 67 of file KalmanTrack_module.cc.

Referenced by CellToPlaneSort(), CellToPlaneSortHits(), CreateHitMaps(), PlaneToCellSort(), PlaneToCellSortHits(), and SegmentFinder().

std::vector<trk::LocatedChan> trk::KalmanTrack::fSliceChans
private

Definition at line 65 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), CreateHitMaps(), FilterOnly(), FilterTracker(), FindTracks(), IsSinglePlane(), produce(), RemoveHitsFromSignal(), SegmentFinder(), and SingleSegment().

std::vector<trk::LocatedChan> trk::KalmanTrack::fSliceChansCellSort
private

Definition at line 66 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), CreateHitMaps(), FilterOnly(), FilterTracker(), FindTracks(), and SingleSegment().

double trk::KalmanTrack::fSysVariance = 0.0001
private

Definition at line 62 of file KalmanTrack_module.cc.

Referenced by FilterOnly(), and FilterTracker().

bool trk::KalmanTrack::fTrySingle
private

Definition at line 71 of file KalmanTrack_module.cc.

Referenced by FindTracks(), and SegmentFinder().

geo::View_t trk::KalmanTrack::fView
private

Definition at line 68 of file KalmanTrack_module.cc.

Referenced by CheckTrack(), FilterOnly(), FilterTracker(), FindTracks(), produce(), SegmentFinder(), and SingleSegment().

int trk::KalmanTrack::fWhoaNelly
private

Definition at line 43 of file KalmanTrack_module.cc.

Referenced by reconfigure(), and SegmentFinder().

int trk::KalmanTrack::NCells[2]
private

Definition at line 63 of file KalmanTrack_module.cc.

Referenced by FilterTracker(), and produce().

The documentation for this class was generated from the following file:
Generated on Wed Apr 21 2021 14:19:18 for NOvA by   doxygen 1.8.11