17 #include "cetlib_except/exception.h" 26 #include "DAQChannelMap/DAQChannelMap.h" 46 const std::vector<double>& zx,
47 const std::vector<double>&
y,
48 const std::vector<double>& zy,
49 std::vector<double>&
w);
50 virtual void FitView(
const std::vector<double>& x,
51 const std::vector<double>& zx,
52 std::vector<double>& wx,
55 const std::vector<double>& tdc,
59 virtual double DsqrToLine(
double x0,
double y0,
61 double x2,
double y2);
63 const std::vector<double>& y,
64 const std::vector<double>& w,
65 double* x1,
double* y1,
66 double* x2,
double* y2);
67 virtual double LinFit(
const std::vector<double>& x,
68 const std::vector<double>& y,
69 const std::vector<double>& w,
70 double* x1,
double* y1,
71 double* x2,
double* y2);
72 virtual double UtilLinFit(
const std::vector<double>& x,
73 const std::vector<double>& y,
74 const std::vector<double>& w,
75 double&
m,
double&
c);
82 std::vector<double>
fR;
83 std::vector<double>
fT;
100 fR[0] = 0.0;
fT[0] = 0.0;
110 produces< std::vector<novaddt::Track> >();
111 produces< std::vector<novaddt::HitList> >();
114 produces< art::Assns<novaddt::Track, novaddt::HitList> >();
134 for(
unsigned int i = 0;
i < slices->size(); ++
i){
140 std::unique_ptr< std::vector<novaddt::Track> > trackcol (
new std::vector<novaddt::Track> );
141 std::unique_ptr< std::vector<novaddt::HitList> > hitListcol(
new std::vector<novaddt::HitList>);
146 for(
size_t i = 0;
i < slices->size(); ++
i){
150 std::vector<double> xzViewX;
151 std::vector<double> xzViewZ;
152 std::vector<double> yzViewY;
153 std::vector<double> yzViewZ;
154 std::vector<double> xzViewW;
155 std::vector<double> yzViewW;
157 std::vector<double> xzTDC;
158 std::vector<double> yzTDC;
162 double lTDC = slicedhits->at(0).TDC().val;
167 for(
auto const&
j : *slicedhits){
172 LOG_DEBUG(
"Tracking") <<
"\thit[" << counter
173 <<
"]: TDC: " <<
j.TDC().val
174 <<
", ADC: " <<
j.ADC().val
175 <<
", plane: " <<
j.Plane().val
177 <<
"-cell: " <<
j.Cell().val
187 xzViewX.push_back(
j.Cell().val);
188 xzViewZ.push_back(
j.Plane().val);
189 xzViewW.push_back(0);
190 xzTDC.push_back(
j.TDC().val);
193 yzViewY.push_back(
j.Cell().val);
194 yzViewZ.push_back(
j.Plane().val);
195 yzViewW.push_back(0);
196 yzTDC.push_back(
j.TDC().val);
205 this->
SeedWeights(xzViewX, xzViewZ, yzViewY, yzViewZ, xzViewW);
206 this->
SeedWeights(yzViewY, yzViewZ, xzViewX, xzViewZ, yzViewW);
209 double xzStart[2] = {0.};
210 double xzEnd[2] = {0.};
211 double yzStart[2] = {0.};
212 double yzEnd[2] = {0.};
214 this->
FitView(xzViewX, xzViewZ, xzViewW, xzStart, xzEnd, xzTDC, eTDC, lTDC);
215 this->
FitView(yzViewY, yzViewZ, yzViewW, yzStart, yzEnd, yzTDC, eTDC, lTDC);
218 if(xzStart[0] == 0. &&
225 yzEnd[1] == 0.)
continue;
227 double xzSlope = (xzStart[1] - xzEnd[1])/(xzStart[0] - xzEnd[0]);
228 double yzSlope = (yzStart[1] - yzEnd[1])/(yzStart[0] - yzEnd[0]);
229 double xzInt = xzStart[1] - xzSlope*xzStart[0];
230 double yzInt = yzStart[1] - yzSlope*yzStart[0];
237 for(
auto const&
h : *slicedhits){
251 if((hl_xz.size() < 2) &&
252 (hl_yz.size() < 2))
continue;
254 LOG_DEBUG(
"Tracking") <<
"--- Track found: x: " << xzStart[1]
256 <<
", z: " << xzStart[0]
258 <<
", n hits: " << hl_xz.size()
260 LOG_DEBUG(
"Tracking") <<
"--- Track found: y: " << yzStart[1]
262 <<
", z: " << yzStart[0]
264 <<
", n hits: " << hl_yz.size()
276 xzStart[1], xzStart[0],startt,
277 xzEnd[1], xzEnd[0],endt
280 yzStart[1], yzStart[0],startt,
281 yzEnd[1], yzEnd[0],endt
284 trackcol->push_back(x_track);
285 hitListcol->push_back(hl_xz);
286 util::CreateAssn(*
this, e, *trackcol, *hitListcol, *assn, hitListcol->
size()-1, hitListcol->size());
288 trackcol->push_back(y_track);
289 hitListcol->push_back(hl_yz);
290 util::CreateAssn(*
this, e, *trackcol, *hitListcol, *assn, hitListcol->
size()-1, hitListcol->size());
296 e.
put(std::move(trackcol));
297 e.
put(std::move(hitListcol));
298 e.
put(std::move(assn));
315 const std::vector<double>& zx,
316 const std::vector<double>&
y,
317 const std::vector<double>& zy,
318 std::vector<double>&
w)
320 for (
unsigned int i=0;
i<zx.size(); ++
i) {
326 for (
unsigned int j=0;
j<zx.size(); ++
j) {
329 d =
sqrt( (zx[
i]-zx[
j])*(zx[
i]-zx[j]) + (x[
i]-x[j])*(x[
i]-x[j]) );
330 if (d<3.0) w[
i] += 1.0;
335 for (
unsigned int j=0;
j<zy.size(); ++
j) {
337 if (d<3.0) w[
i] += 0.1;
346 const std::vector<double>& zx,
347 std::vector<double>& wx,
350 const std::vector<double>& tdc,
364 double intercept = 0;
367 for (
unsigned int i=0;
i<
fR.size(); ++
i) {
368 static const double eps = 1.0E-6;
372 for (
unsigned int j=0;
j<wx.size(); ++
j) {
373 double d = this->
DsqrToLine(zx[
j], x[j], zx1, x1, zx2, x2);
382 if (wx[j]<eps) wx[
j] = eps;
387 if(zx.size() != x.size() ||
388 wx.size() != x.size() ||
396 double dzx = zx2-zx1;
397 if (dzx<=0.0) dzx = 1.0E-9;
400 intercept = (x1*zx2-x2*zx1)/dzx;
413 for (
size_t i = 0;
i < wx.size(); ++
i) {
415 if (zx[
i] < start[0]) start[0] = zx[
i];
416 if (zx[
i] > end[0]) end[0] = zx[
i];
417 if (tdc[
i] < eTDC) eTDC = tdc[
i];
418 if (tdc[
i] > lTDC) lTDC = tdc[
i];
422 LOG_DEBUG(
"Tracking") <<
"start: " << start[0] <<
", end: " << end[0] <<
", slope: " << slope <<
", intercept: " << intercept <<
std::endl;
423 start[1] = slope*start[0] + intercept;
424 end[1] = slope*end[0] + intercept;
427 if (start[1]<=0.) start[1] = 1.0E-9;
428 if (end[1]<=0.) end[1] = 1.0E-9;
429 if (start[1]>383.) start[1] = 383.;
430 if (end[1]>383.) end[1] = 383.;
437 double x1,
double y1,
438 double x2,
double y2)
442 if (x1==x2 && y1==y2)
return (x0-x1)*(x0-
x1) +(y0-y1)*(y0-
y1);
445 ((x2-x1)*(y1-y0)-(x1-x0)*(y2-
y1))*((x2-
x1)*(y1-y0)-(x1-x0)*(y2-y1)) /
446 ( (x2-x1)*(x2-
x1) + (y2-y1)*(y2-
y1) );
452 const std::vector<double>&
y,
453 const std::vector<double>&
w,
454 double*
x1,
double*
y1,
455 double*
x2,
double*
y2)
461 unsigned register int i;
464 double xlo = DBL_MAX;
465 double xhi = -DBL_MAX;
466 double ylo = DBL_MAX;
467 double yhi = -DBL_MAX;
468 for (i=0; i<w.size(); ++
i) {
469 if(w[i] == 0)
continue;
470 if (x[i]<xlo) xlo = x[
i];
471 if (y[i]<ylo) ylo = y[
i];
472 if (x[i]>xhi) xhi = x[
i];
473 if (y[i]>yhi) yhi = y[
i];
490 double xx1, yy1, xx2, yy2;
508 for (i=0; i<w.size(); ++
i) {
512 Swxy += w[
i]*x[
i]*y[
i];
513 Swy2 += w[
i]*y[
i]*y[
i];
514 Swx2 += w[
i]*x[
i]*x[
i];
521 const double SwSq = Sw*Sw;
522 const double SwxSq = Swx*Swx;
523 const double SwySq = Swy*Swy;
539 if (C>0.0) C =
sqrt(C);
541 double D = -Swx*Swy+Swxy*Sw;
555 M1 = -(-SwxSq + Swx2*Sw + SwySq - Swy2*Sw -
C)/(2.0*D);
556 B1 = -(-Swy+M1*Swx)/Sw;
570 for (i=0; i<w.size(); ++
i) {
571 chi2 += w[
i]*this->
DsqrToLine(x[i],y[i],*x1,*y1,*x2, *y2);
575 double chi2lin, linx1, linx2, liny1, liny2;
577 chi2lin =
LinFit(x,y,w,&linx1,&liny1,&linx2,&liny2);
596 const std::vector<double>&
y,
597 const std::vector<double>&
w,
598 double*
x1,
double*
y1,
599 double*
x2,
double*
y2)
601 register unsigned int i;
604 assert(x.size()==y.size());
605 assert(x.size()==w.size());
614 for (i=1; i<w.size(); ++
i) {
615 if (x[i]<xlo) xlo = x[
i];
616 if (x[i]>xhi) xhi = x[
i];
617 if (y[i]<ylo) ylo = y[
i];
618 if (y[i]>yhi) yhi = y[
i];
621 if (yhi-ylo==0.0 && xhi-xlo==0.0) {
623 <<
" All points identical in line fit!" <<
std::endl;
635 if (yhi-ylo > xhi-xlo) {
636 return this->
LinFit(y,x,w,y1,x1,y2,x2);
653 const std::vector<double>&
y,
654 const std::vector<double>&
w,
655 double&
m,
double&
c)
658 assert(x.size() == y.size());
659 assert(x.size() == w.size());
669 for(
unsigned int i = 0;
i < w.size(); ++
i) {
673 Swx2 += w[
i]*x[
i]*x[
i];
674 Swxy += w[
i]*x[
i]*y[
i];
675 Swy2 += w[
i]*y[
i]*y[
i];
677 const double d = Sw*Swx2 - Swx*Swx;
678 m = (Sw*Swxy - Swx*Swy)/d;
679 c = (Swy*Swx2 - Swx*Swxy)/d;
682 Swy2 - 2.0*m*Swxy - 2.0*c*Swy + 2.0*m*c*Swx +
#define LOG_DEBUG(stream)
virtual void produce(art::Event &e)
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.
Float_t y1[n_points_granero]
virtual double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Float_t x1[n_points_granero]
std::vector< DAQHit > HitList
::xsd::cxx::tree::exception< char > exception
virtual void FitView(const std::vector< double > &x, const std::vector< double > &zx, std::vector< double > &wx, double *start, double *End, const std::vector< double > &tdc, double &eTDC, double &lTDC)
DEFINE_ART_MODULE(TestTMapFile)
std::vector< double > fT
Transition width of weight function.
virtual void SeedWeights(const std::vector< double > &x, const std::vector< double > &zx, const std::vector< double > &y, const std::vector< double > &zy, std::vector< double > &w)
virtual double UtilLinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Identifier for the Y measuring view of the detector (side)
ProductID put(std::unique_ptr< PROD > &&product)
void push_back(Ptr< U > const &p)
Identifier for the X measuring view of the detector (top)
fvar< T > exp(const fvar< T > &x)
std::string fSliceInstance
virtual double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
EventNumber_t event() const
std::string _assn_to_slice_instance
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
std::vector< double > fR
Range of weight function.
TrackFit(fhicl::ParameterSet const &p)
assert(nhit_max >=nhit_nbins)
virtual 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)