Public Types | Public Member Functions | Static 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::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

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

Public Member Functions

 KalmanTrack (fhicl::ParameterSet const &pset)
 
virtual ~KalmanTrack ()
 
void produce (art::Event &evt)
 
void reconfigure (const fhicl::ParameterSet &p)
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< T > consumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< T > consumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< T > mayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< T > mayConsumeView (InputTag const &it)
 
base_engine_t & createEngine (seed_t seed)
 
base_engine_t & createEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_t & createEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumer > non_module_context ()
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

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 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

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

Definition at line 35 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  {
149  reconfigure(pset);
150  produces< std::vector<rb::Track> >();
151  produces< art::Assns<rb::Track, rb::Cluster> >();
152  }
trk::KalmanTrack::reconfigure
void reconfigure(const fhicl::ParameterSet &p)
Definition: KalmanTrack_module.cc:161
trk::KalmanTrack::~KalmanTrack ( )
virtual

Definition at line 155 of file KalmanTrack_module.cc.

156  {
157  }

Member Function Documentation

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

Definition at line 2448 of file KalmanTrack_module.cc.

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

Referenced by FindTracks().

2449  {
2450  std::vector<std::array<double, 2>> planeSortedVec(cellSortedVec.size());
2451  // loop over the cellSortedVec and set the corresponding planeSortedHit to true
2452  for(unsigned int i = 0; i < cellSortedVec.size(); ++i){
2453  int planesort = 0;
2454  // get the corresponding plane sorted hit location
2455  for(unsigned int j = 0; j < fPlaneToCell.size(); ++j){
2456  if(fPlaneToCell[j] == (int)i){
2457  planesort = j;
2458  break;
2459  }
2460  }
2461  planeSortedVec[planesort] = cellSortedVec[i];
2462  }
2463  cellSortedVec = planeSortedVec;
2464  }
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 2416 of file KalmanTrack_module.cc.

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

Referenced by FindTracks().

2417  {
2418  std::vector<bool> planeSortedHits(cellSortedHits.size());
2419  // loop over the cellSortedHits and set the corresponding planeSortedHit to true
2420  for(unsigned int i = 0; i < cellSortedHits.size(); ++i){
2421  int planesort = 0;
2422  // get the corresponding plane sorted hit location
2423  for(unsigned int j = 0; j < fPlaneToCell.size(); ++j){
2424  if(fPlaneToCell[j] == (int)i){
2425  planesort = j;
2426  break;
2427  }
2428  }
2429  planeSortedHits[planesort] = cellSortedHits[i];
2430  }
2431  cellSortedHits = planeSortedHits;
2432  }
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 2115 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().

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

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
it
set< int >::iterator it
Definition: ProduceBadRunsByFailureMode.C:101
art::ProductToken::invalid
static ProductToken< T > invalid()
Definition: ProductToken.h:47
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:137
T
double T
Definition: Xdiff_gwt.C:5
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:137
T
double T
Definition: Xdiff_gwt.C:5
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
it
set< int >::iterator it
Definition: ProduceBadRunsByFailureMode.C:101
art::ViewToken::invalid
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:137
T
double T
Definition: Xdiff_gwt.C:5
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:135
int trk::KalmanTrack::CountHits ( const std::vector< bool > &  trkHits)
private

Definition at line 2467 of file KalmanTrack_module.cc.

References fNChans, and MECModelEnuComparisons::i.

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

2468  {
2469  int nhits = 0;
2470  for(int i = 0; i < fNChans; ++i){
2471  if(trkHits[i]){ ++nhits; }
2472  }
2473  return nhits;
2474  }
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited

Referenced by evgen::CosmicPionGen::CosmicPionGen(), evgen::CosmicsGen::CosmicsGen(), fuzz::FuzzyKVertex::FuzzyKVertex(), g4n::G4Gen::G4Gen(), photrans::ImprovedTransport::ImprovedTransport(), art::MixFilter< T >::initEngine_(), evgen::IsoGen::IsoGen(), rsim::MakeNoiseSpectrumFile::MakeNoiseSpectrumFile(), photrans::PhotonSystRewriter::PhotonSystRewriter(), photrans::PhotonTransport::PhotonTransport(), rsim::ReadoutSim::ReadoutSim(), evgen::SingleGen::SingleGen(), and evgen::SupernovaGen::SupernovaGen().

base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
void trk::KalmanTrack::CreateHitMaps ( )
private

Definition at line 2386 of file KalmanTrack_module.cc.

References fNChans, fPlaneToCell, fSliceChans, and fSliceChansCellSort.

Referenced by FindTracks().

2387  {
2388  fPlaneToCell.resize(fNChans);
2389  // loop over all the plane sorted hits
2390  for(int ipl = 0; ipl < fNChans; ++ipl){
2391  art::Ptr<rb::CellHit> plhit = fSliceChans[ipl].GetHit();
2392  // loop over all the cell sorted hits and compare to this hit
2393  for(unsigned int icl = 0; icl < fSliceChansCellSort.size(); ++icl){
2394  if(plhit == fSliceChansCellSort[icl].GetHit()){
2395  fPlaneToCell[ipl] = icl;
2396  break;
2397  }
2398  }
2399  }
2400  } // 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
CurrentProcessingContext const* art::EDProducer::currentContext ( ) const
protectedinherited
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 1754 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().

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

Referenced by FindTracks().

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

Definition at line 398 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(), getGoodRuns4SAM::n, 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().

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

Referenced by art::MixFilter< T >::initEngine_().

int trk::KalmanTrack::GetFirstHit ( const std::vector< bool > &  trkHits)
private

Definition at line 2502 of file KalmanTrack_module.cc.

References fNChans, and MECModelEnuComparisons::i.

Referenced by FilterOnly(), and FindTracks().

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

Definition at line 2515 of file KalmanTrack_module.cc.

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

2516  {
2517  int lastHit = -1;
2518  for(int i = fNChans-1; i >= 0; --i){
2519  if(trkHits[i]){
2520  lastHit = i;
2521  break;
2522  }
2523  }
2524  return lastHit;
2525  }
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

Referenced by skim::NueSkimmer::CopyMichelSlice(), and skim::NueSkimmer::CopyMichelTrack().

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
art::EDProducer::moduleDescription_
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
fetch_tb_beamline_files.md
md
Definition: fetch_tb_beamline_files.py:11
bool trk::KalmanTrack::IsSinglePlane ( const std::vector< bool > &  trkHits)
private

Definition at line 2477 of file KalmanTrack_module.cc.

References fNChans, fSliceChans, and MECModelEnuComparisons::i.

Referenced by FindTracks(), and SingleSegment().

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

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
it
set< int >::iterator it
Definition: ProduceBadRunsByFailureMode.C:101
art::ProductToken::invalid
static ProductToken< T > invalid()
Definition: ProductToken.h:47
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:137
T
double T
Definition: Xdiff_gwt.C:5
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:137
T
double T
Definition: Xdiff_gwt.C:5
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
it
set< int >::iterator it
Definition: ProduceBadRunsByFailureMode.C:101
art::ViewToken::invalid
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
art::Consumer::consumables_
ConsumableProducts consumables_
Definition: Consumer.h:137
T
double T
Definition: Xdiff_gwt.C:5
art::Consumer::moduleContext_
bool moduleContext_
Definition: Consumer.h:135
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID(), and string.

41  {
42  return true;
43  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited

Referenced by art::RootOutput::beginSubRun(), art::ProcessPackage< L >::postScheduleSignal(), art::BeginEndPackage< Level::Run >::Begin::postScheduleSignal(), art::BeginEndPackage< Level::Run >::End::postScheduleSignal(), art::BeginEndPackage< Level::SubRun >::Begin::postScheduleSignal(), art::BeginEndPackage< Level::SubRun >::End::postScheduleSignal(), art::ProcessPackage< L >::preScheduleSignal(), art::BeginEndPackage< Level::Run >::Begin::preScheduleSignal(), art::BeginEndPackage< Level::SubRun >::Begin::preScheduleSignal(), art::EmptyEvent::readRun_(), and art::EmptyEvent::readSubRun_().

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

Definition at line 2435 of file KalmanTrack_module.cc.

References fPlaneToCell, and MECModelEnuComparisons::i.

Referenced by FindTracks().

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

Definition at line 2403 of file KalmanTrack_module.cc.

References fPlaneToCell, and MECModelEnuComparisons::i.

Referenced by FindTracks(), and SingleSegment().

2404  {
2405  std::vector<bool> cellSortedHits(planeSortedHits.size());
2406  // loop over the planeSortedHits and set the corresponding cellSortedHit to true
2407  for(unsigned int i = 0; i < planeSortedHits.size(); ++i){
2408  // get the corresponding cell sorted hit location
2409  int cellsort = fPlaneToCell[i];
2410  cellSortedHits[cellsort] = planeSortedHits[i];
2411  }
2412  planeSortedHits = cellSortedHits;
2413  }
trk::KalmanTrack::fPlaneToCell
std::vector< int > fPlaneToCell
Definition: KalmanTrack_module.cc:67
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
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 181 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(), stan::math::floor(), 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::Event::put(), art::PtrVector< T >::reserve(), trk::LocatedChan::SetHit(), art::PtrVector< T >::size(), trk::SortByPlane(), and geo::PlaneGeo::View().

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

Definition at line 161 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().

162  {
163  fMaxSliceHits = pset.get< double >("MaxSliceHits") ;
164  fMaxHitCut = pset.get< double >("MaxHitCut") ;
165  fWhoaNelly = pset.get< int >("WhoaNelly") ;
166  fClusterInput = pset.get< std::string >("ClusterInput") ;
167  fDeltaChiAllowed = pset.get< double >("DeltaChiAllowed") ;
168  fMinHits = pset.get< int >("MinHits") ;
169  fGapAllowed = pset.get< int >("GapAllowed") ;
170  fLargeSliceHits = pset.get< int >("LargeSliceHits") ;
171  fBadChannels = pset.get< bool >("BadChannels") ;
172  fLongTrackLoop = pset.get< bool >("LongTrackLoop") ;
173  fMinGapProb = pset.get< double >("MinGapProb") ;
174  fObeyPreselection = pset.get< bool >("ObeyPreselection");
175  }
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 1726 of file KalmanTrack_module.cc.

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

Referenced by FindTracks().

1727  {
1728 
1729  /// Following is added for new hit storage method
1730  std::vector<trk::LocatedChan> newfSignalChan;
1731  // loop over all the signal hits
1732  for(int i = 0; i<fNChans; ++i){
1733  // get current hit
1734  art::Ptr<rb::CellHit> chanhit = fSliceChans[i].GetHit();
1735  // define a boolean to track if signal hit is in track hit
1736  bool inTrack = false;
1737  for(unsigned int j = 0; j < track.NCell(); ++j){
1738  // get track hit
1739  art::Ptr<rb::CellHit> thit = track.Cell(j);
1740  if(chanhit == thit){
1741  inTrack = true;
1742  break;
1743  }
1744  }
1745  if(!inTrack){ newfSignalChan.push_back(fSliceChans[i]); }
1746  }
1747  fNChans = newfSignalChan.size();
1748  fSliceChans.resize(fNChans);
1749  fSliceChans = newfSignalChan;
1750 
1751  }
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 1666 of file KalmanTrack_module.cc.

References fNChans.

Referenced by FindTracks().

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

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

Referenced by art::RootOutput::endJob().

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 996 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().

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

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 Sat Nov 28 2020 05:22:15 for NOvA by   doxygen 1.8.11