14 #include "Utilities/func/MathUtil.h" 36 std::vector<rb::Track> tracks;
50 std::vector<double> xzViewX;
51 std::vector<double> xzViewZ;
52 std::vector<double> yzViewY;
53 std::vector<double> yzViewZ;
54 std::vector<double> xzViewW;
55 std::vector<double> yzViewW;
57 const int hitMax = slice->
NCell();
58 for(
int hitIdx = 0; hitIdx < hitMax; ++hitIdx){
64 xzViewX.push_back(xyz[0]);
65 xzViewZ.push_back(xyz[2]);
66 xzViewW.push_back(sig);
69 <<
"x view " << xzViewZ[xzViewZ.size()-1]
70 <<
" " << xzViewX[xzViewX.size()-1]
71 <<
" " << chit->
TNS();
74 yzViewY.push_back(xyz[1]);
75 yzViewZ.push_back(xyz[2]);
76 yzViewW.push_back(sig);
79 <<
"y view " << yzViewZ[yzViewZ.size()-1]
80 <<
" " << yzViewY[yzViewY.size()-1]
81 <<
" " << chit->
TNS();
85 if(xzViewZ.size() < 2 && yzViewZ.size() < 2){
87 <<
" Not enough hits in either view to fit for track";
97 double xzStart[2] = {0.};
98 double xzEnd[2] = {0.};
99 unsigned int xzNDrop = this->
FitView(xzViewZ, xzViewX, xzViewW, xzChi2, xzStart, xzEnd);
100 if(xzEnd[0]==xzStart[0] && xzEnd[1]==xzStart[1] && xzViewW.size()>1) flagxz =
false ;
101 if(xzViewZ.size() == 1){
102 xzStart[0] = xzEnd[0] = xzViewZ[0];
103 xzStart[1] = xzEnd[1] = xzViewX[0];
107 double yzStart[2] = {0.};
108 double yzEnd[2] = {0.};
109 unsigned int yzNDrop = this->
FitView(yzViewZ, yzViewY, yzViewW, yzChi2, yzStart, yzEnd);
110 if(yzEnd[0]==yzStart[0] && yzEnd[1]==yzStart[1] && yzViewW.size()>1) flagyz =
false ;
111 if(yzViewZ.size() == 1){
112 yzStart[0] = yzEnd[0] = yzViewZ[0];
113 yzStart[1] = yzEnd[1] = yzViewY[0];
117 <<
"# xz hits = " << xzViewZ.size()-xzNDrop
119 <<
" # yz hits = " << yzViewZ.size()-yzNDrop
120 <<
" chi^2 = " << yzChi2;
128 if(
fabs(xzEnd[0]-xzStart[0]) < 1){
130 double ZaveY = yzEnd[0]-yzStart[0]/2.;
132 if(
fabs(xzStart[1]) <
fabs(xzEnd[1]) ) {
133 if( xzStart[0] < ZaveY ) {
136 if( xzStart[0] > ZaveY ) {
141 if( xzEnd[0] > ZaveY ) {
144 if( xzEnd[0] < ZaveY ) {
150 if(
fabs(yzEnd[0]-yzStart[0]) < 1){
152 double ZaveX = xzEnd[0]-xzStart[0]/2.;
154 if(
fabs(yzStart[1]) <
fabs(yzEnd[1]) ) {
155 if( yzStart[0] < ZaveX ) {
158 if( yzStart[0] > ZaveX ) {
163 if( yzEnd[0] > ZaveX ) {
166 if( yzEnd[0] < ZaveX ) {
173 double startz = TMath::Min(yzStart[0], xzStart[0]);
174 double endz = TMath::Max(yzEnd[0], xzEnd[0]);
178 if(xzViewZ.size() < 2){
183 if(yzViewZ.size() < 2){
193 const double xzSlope = (xzEnd[1]-xzStart[1])/(xzEnd[0]-xzStart[0]);
194 xzStart[1] += (startz-xzStart[0])*xzSlope;
195 xzEnd[1] += (endz-xzEnd[0])*xzSlope;
197 const double yzSlope = (yzEnd[1]-yzStart[1])/(yzEnd[0]-yzStart[0]);
198 yzStart[1] += (startz-yzStart[0])*yzSlope;
199 yzEnd[1] += (endz-yzEnd[0])*yzSlope;
202 xzStart[0] = yzStart[0] = startz;
203 xzEnd[0] = yzEnd[0] = endz;
207 << xzStart[0] <<
" " << xzStart[1] <<
" " 208 << xzEnd[0] <<
" " << xzEnd[1] <<
std::endl 209 << yzStart[0] <<
" " << yzStart[1] <<
" " 210 << yzEnd[0] <<
" " << yzEnd[1] <<
std::endl;
230 for(
unsigned int hitIdx = 0; hitIdx < slice->
NXCell(); ++hitIdx){
231 if(xzViewW[hitIdx] != 0){
234 if (ch->
Plane() < minPlane) minPlane = ch->
Plane();
238 for(
unsigned int hitIdx = 0; hitIdx < slice->
NYCell(); ++hitIdx){
239 if(yzViewW[hitIdx] != 0){
242 if (ch->
Plane() < minPlane) minPlane = ch->
Plane();
249 <<
"No hits left on the track after clean up stage ";
258 TVector3
start(xzStart[1], yzStart[1], startz);
259 TVector3 stop(xzEnd[1], yzEnd[1], endz);
261 if(start[1] < stop[1]) {
278 double minZ =
std::min(start[2],stop[2]);
279 double maxZ =
std::max(start[2],stop[2]);
287 int prevPlane = chits.
front()->Plane();
288 double xyz[3] = {0.};
289 for(
size_t c = 1;
c < chits.
size()-1; ++
c){
293 if ((plane == minPlane) || (plane == maxPlane))
continue;
295 if(plane != prevPlane) {
296 geom->
Plane(plane)->
Cell(chits.
at(
c)->Cell())->GetCenter(xyz);
298 if ((xyz[2] > minZ) && (xyz[2] < maxZ)){
299 TVector3
dir = track.
Dir().Unit();
300 TVector3 point(start.X() + (xyz[2] - start.Z())*dir.X()/dir.Z(),
301 start.Y() + (xyz[2] - start.Z())*dir.Y()/dir.Z(),
312 tracks.push_back(track);
320 if(flagxz&&!flagyz) {
324 for(
unsigned int hitIdx = 0; hitIdx < slice->
NXCell(); ++hitIdx){
330 tracks.push_back(track);
334 if(flagyz&&!flagxz) {
337 for(
unsigned int hitIdx = 0; hitIdx < slice->
NYCell(); ++hitIdx){
343 tracks.push_back(track);
352 std::vector<double> &
y,
353 std::vector<double> &w_orig,
359 std::vector<double>
w = w_orig;
361 const unsigned int nhit = w.size();
370 unsigned int ndrop = 0;
371 for(
unsigned int i = 0;
i <
nhit; ++
i){
372 double closest = 1e10;
373 for(
unsigned int j = 0;
j <
nhit; ++
j){
386 if(nhit-ndrop == 0)
return 0;
391 start[1] = end[1] = y[0];
401 if(start[0] > end[0]){
408 bool droppedAny =
true;
411 &start[0], &start[1],
414 LOG_DEBUG(
"CosmicTrack") <<
"end points from fit " 415 << start[0] <<
" " << start[1] <<
" " << end[0] <<
" " << end[1];
424 for(
unsigned int n = 0;
n <
nhit; ++
n){
426 if(w[
n] == 0)
continue;
439 LOG_DEBUG(
"CosmicTrack") <<
"dropping hit at (" 440 << x[dSqMaxIdx] <<
", " << y[dSqMaxIdx] <<
")" 441 <<
" dSq = " << dSqMax;
450 bool addedAny =
false;
452 for(
unsigned int n = 0;
n <
nhit; ++
n){
454 if(w[
n] != 0)
continue;
461 <<
"adding back hit at (" << x[
n] <<
", " << y[
n]
462 <<
")" <<
" dSq = " << dSq;
473 &start[0], &start[1],
477 LOG_DEBUG(
"CosmicTrack") <<
"end points from fit " 478 << start[0] <<
" " << start[1] <<
" " 479 << end[0] <<
" " << end[1];
484 if(
fabs(start[0]-end[0])>0.01&&
fabs(start[1]-end[1])>0.01) {
486 const TVector3 slope(end[0]-start[0],end[1]-start[1],0);
487 const TVector3 intercept(0,(-end[0]/slope[0])*slope[1]+end[1],0);
488 double tmin[2], tmax[2];
495 for(
unsigned int n = 0;
n <
nhit; ++
n){
496 if(w[
n] == 0)
continue;
498 TVector3 point(x[
n],y[n],0);
503 if(closest[0]<tmin[0]) {
504 tmin[0] = closest[0];
505 tmin[1] = closest[1];
507 if(closest[0]>tmax[0]) {
508 tmax[0] = closest[0];
509 tmax[1] = closest[1];
513 if(
fabs(tmin[0]-end[0])<
fabs(tmin[0]-start[0])) {
538 std::vector<double> zs, ts, ws;
540 const int hitMax = slice->
NCell();
541 for(
int hitIdx = 0; hitIdx < hitMax; ++hitIdx){
546 zs.push_back(xyz[2]);
547 ts.push_back(chit->
TNS());
552 double ztStart[2] = {0.};
553 double ztEnd[2] = {0.};
555 &ztStart[0], &ztStart[1],
556 &ztEnd[0], &ztEnd[1]);
560 const double ztSlope = (ztEnd[1]-ztStart[1])/(ztEnd[0]-ztStart[0]);
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
An algorithm to perform cosmic ray track fitting.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
fvar< T > fabs(const fvar< T > &x)
unsigned short Plane() const
const CellGeo * Cell(int icell) const
Vertical planes which measure X.
double fDHitGood
Maximum distance from the line for a hit to be considered part of it.
A collection of associated CellHits.
A rb::Prong with full reconstructed trajectory.
bool ClampRayToDetectorVolume(TVector3 *p0, TVector3 *p1, const GeometryBase *geom)
If either endpoint is outside the detector move it to the edge.
const PlaneGeo * Plane(unsigned int i) const
Track finder for cosmic rays.
CosmicTrackAlg(fhicl::ParameterSet const &pset)
T sqr(T x)
More efficient square function than pow(x,2)
Horizontal planes which measure Y.
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
unsigned short Cell() const
Track finder for cosmic rays.
virtual ~CosmicTrackAlg()
void push_back(Ptr< U > const &p)
double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
In two dimensions, return the perpendicular distance from a point (x0,y0) to a line defined by end po...
void AppendTrajectoryPoint(TVector3 pt)
T get(std::string const &key) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
Collect Geo headers and supply basic geometry functions.
reference at(size_type n)
unsigned int FitView(std::vector< double > &x, std::vector< double > &y, std::vector< double > &w, double &chiSqr, double *start, double *end)
unsigned int NYCell() const
Number of cells in the y-view.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
unsigned int NXCell() const
Number of cells in the x-view.
assert(nhit_max >=nhit_nbins)
bool IsTrackDownstreamFromTiming(const rb::Cluster *slice) const
T min(const caf::Proxy< T > &a, T b)
Encapsulate the geometry of one entire detector (near, far, ndos)
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
std::vector< rb::Track > MakeTrack(const rb::Cluster *slice)