32 #include "Utilities/AssociationUtil.h" 55 const std::vector<double>&
y,
56 const std::vector<double>&
w,
57 double& sigXY,
double& sigZ,
58 double&
x1,
double&
y1,
59 double&
x2,
double&
y2,
60 double& sigy1,
double& sigy2 );
155 fGeoTree = tfs->
make<TTree>(
"GeometryInfo",
"Geometry Info");
156 fOutTree = tfs->
make<TTree>(
"SplitTracks",
"Split Tracks");
166 fOutTree->Branch(
"trk_ID", &
trk_ID);
170 fOutTree->Branch(
"trk_StopX", &
trk_StopX);
171 fOutTree->Branch(
"trk_StopY", &
trk_StopY);
172 fOutTree->Branch(
"trk_StopZ", &
trk_StopZ);
188 subtrkXYstart = tfs->
make<TH2D>(
"subtrkXYstart",
"Downstream Subtrack Start (XY);X-view Cell;Y-view Cell",100,-775,775,100,-775,775);
189 subtrkXZstart = tfs->
make<TH2D>(
"subtrkXZstart",
"Downstream Subtrack Start (XZ);Plane;X-view Cell",600,0,6000,1550,-775,775);
190 subtrkYZstart = tfs->
make<TH2D>(
"subtrkYZstart",
"Downstream Subtrack Start (YZ);Plane;Y-view Cell",600,0,6000,1550,-775,775);
192 xzViewTrkRes = tfs->
make<TH1D>(
"xzViewTrkRes",
"Track Residual (X-View);Residual",100,-10,10);
193 yzViewTrkRes = tfs->
make<TH1D>(
"yzViewTrkRes",
"Track Residual (Y-View);Residual",100,-10,10);
194 xzViewTrkResVsX = tfs->
make<TH2D>(
"xzViewTrkResVsX",
"Track Residual (X-View) vs Cell Position (X);Cell Number;Residual",32*12,0,32*12,100,-10,10);
195 xzViewTrkResVsZ = tfs->
make<TH2D>(
"xzViewTrkResVsZ",
"Track Residual (X-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
196 yzViewTrkResVsY = tfs->
make<TH2D>(
"yzViewTrkResVsY",
"Track Residual (Y-View) vs Cell Position (Y);Cell Number;Residual",32*12,0,32*12,100,-10,10);
197 yzViewTrkResVsZ = tfs->
make<TH2D>(
"yzViewTrkResVsZ",
"Track Residual (Y-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
199 xzViewSubtrkRes = tfs->
make<TH1D>(
"xzViewSubtrkRes",
"SubTrack Residual (X-View);Residual",100,-10,10);
200 yzViewSubtrkRes = tfs->
make<TH1D>(
"yzViewSubtrkRes",
"SubTrack Residual (Y-View);Residual",100,-10,10);
201 xzViewSubtrkResVsX = tfs->
make<TH2D>(
"xzViewSubtrkResVsX",
"SubTrack Residual (X-View) vs Cell Position (X);Cell Number;Residual",32*12,0,32*12,100,-10,10);
202 xzViewSubtrkResVsZ = tfs->
make<TH2D>(
"xzViewSubtrkResVsZ",
"SubTrack Residual (X-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
203 yzViewSubtrkResVsY = tfs->
make<TH2D>(
"yzViewSubtrkResVsY",
"SubTrack Residual (Y-View) vs Cell Position (Y);Cell Number;Residual",32*12,0,32*12,100,-10,10);
204 yzViewSubtrkResVsZ = tfs->
make<TH2D>(
"yzViewSubtrkResVsZ",
"SubTrack Residual (Y-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
206 xzViewRefitRes = tfs->
make<TH1D>(
"xzViewRefitRes",
"Refitted SubTrack Residual (X-View);Residual",100,-10,10);
207 yzViewRefitRes = tfs->
make<TH1D>(
"yzViewRefitRes",
"Refitted SubTrack Residual (Y-View);Residual",100,-10,10);
208 xzViewRefitResVsX = tfs->
make<TH2D>(
"xzViewRefitResVsX",
"Refitted SubTrack Residual (X-View) vs Cell Position (X);Cell Number;Residual",32*12,0,32*12,100,-10,10);
209 xzViewRefitResVsZ = tfs->
make<TH2D>(
"xzViewRefitResVsZ",
"Refitted SubTrack Residual (X-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
210 yzViewRefitResVsY = tfs->
make<TH2D>(
"yzViewRefitResVsY",
"Refitted SubTrack Residual (Y-View) vs Cell Position (Y);Cell Number;Residual",32*12,0,32*12,100,-10,10);
211 yzViewRefitResVsZ = tfs->
make<TH2D>(
"yzViewRefitResVsZ",
"Refitted SubTrack Residual (Y-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
215 refitDCAtoCell = tfs->
make<TH1D>(
"refitDCAtoCell",
"Refitted SubTrack DCA to Cell",100,0,10);
253 if ( CosmicTracks->empty() ) {
259 std::unique_ptr< std::vector<rb::Track> > subTrackColl(
new std::vector<rb::Track>);
260 for (
size_t trkItr = 0; trkItr < CosmicTracks->size(); ++trkItr ) {
263 if ( !trackIn->
Is3D() )
continue;
276 if ( maxBlock < minBlock ) {
296 if (!(thruTop&&thruBot) )
continue;
308 std::vector<rb::Cluster> subClusters( maxBlock-minBlock+1,
rb::Cluster() );
310 for (
size_t cellItr = 0; cellItr < trackIn->
NCell(); ++cellItr ) {
314 subClusters.at(block-minBlock).Add(cell);
320 double cellLo[3], cellHi[3];
322 double dOnTrack, dInCell;
323 TVector3 pcaOnTrack, pcaInCell;
324 double dca =
geo::ClosestApproach( trackIn->
Start(), trackIn->
Stop(), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
326 double cellPos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
327 double trkX(0), trkY(0);
330 double res = trkX-cellPos[0];
336 double res = trkY-cellPos[1];
353 for (
size_t clusItr = 0; clusItr < subClusters.size(); ++clusItr ) {
356 std::vector<rb::Track> subTracks =
fTrkAlg->
MakeTrack(&subClusters.at(clusItr));
359 if(subTracks.size() < 1 || subTracks.size() > 1)
continue;
366 if ( subTrack.
NXCell() < 2 )
continue;
367 if ( subTrack.
NYCell() < 2 )
continue;
368 subTrackColl->push_back(subTrack);
371 std::vector<double> xzViewX, xzViewZ, yzViewY, yzViewZ, xzViewW, yzViewW;
372 double zMin = 10000000;
374 for (
size_t i = 0;
i < subTrack.
NCell(); ++
i ) {
381 double cellLo[3], cellHi[3];
383 double dOnTrack, dInCell;
384 TVector3 pcaOnTrack, pcaInCell;
385 double dca =
geo::ClosestApproach( subTrack.
Start(), subTrack.
Stop(), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
387 double pos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
390 double subtrkX(0), subtrkY(0);
393 double res = subtrkX-pos[0];
397 xzViewX.push_back(pos[0]);
398 xzViewZ.push_back(pos[2]);
399 xzViewW.push_back(1.);
402 double res = subtrkY-pos[1];
406 yzViewY.push_back(pos[1]);
407 yzViewZ.push_back(pos[2]);
408 yzViewW.push_back(1.);
410 if ( pos[2] < zMin ) zMin = pos[2];
411 if ( pos[2] > zMax ) zMax = pos[2];
415 double sigX = 3./
sqrt(12);
416 double sigY = 3./
sqrt(12);
417 double sigZ = 5./
sqrt(12);
418 double xStart, yStart, zStart;
419 double xStartErr, yStartErr, zStartErr;
420 double xStop, yStop, zStop;
421 double xStopErr, yStopErr, zStopErr;
422 zStart = zMin; zStartErr = 0.;
423 zStop = zMax; zStopErr = 0.;
424 this->
DemingRegFit( xzViewZ, xzViewX, xzViewW, sigZ, sigX, zStart, xStart, zStop, xStop, xStartErr, xStopErr );
425 this->
DemingRegFit( yzViewZ, yzViewY, yzViewW, sigZ, sigY, zStart, yStart, zStop, yStop, yStartErr, yStopErr );
428 for (
size_t i = 0;
i < subTrack.
NXCell(); ++
i ) {
435 double cellLo[3], cellHi[3];
437 double dOnTrack, dInCell;
438 TVector3 pcaOnTrack, pcaInCell;
439 double dca =
geo::ClosestApproach( TVector3(xStart,yStart,zStart), TVector3(xStop,yStop,zStop), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
441 double pos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
442 double slopeX = (xStop-xStart)/(zStop-zStart);
443 double extrapX = slopeX*(pos[2]-zStart)+xStart;
444 double res = extrapX-pos[0];
449 for (
size_t i = 0;
i < subTrack.
NYCell(); ++
i ) {
456 double cellLo[3], cellHi[3];
458 double dOnTrack, dInCell;
459 TVector3 pcaOnTrack, pcaInCell;
460 double dca =
geo::ClosestApproach( TVector3(xStart,yStart,zStart), TVector3(xStop,yStop,zStop), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
462 double pos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
463 double slopeY = (yStop-yStart)/(zStop-zStart);
464 double extrapY = slopeY*(pos[2]-zStart)+yStart;
465 double res = extrapY-pos[1];
527 const std::vector<double>&
y,
528 const std::vector<double>&
w,
529 double& sigX,
double& sigY,
530 double&
x1,
double&
y1,
531 double&
x2,
double&
y2,
532 double& sigy1,
double& sigy2 )
534 double delta = sigY / sigX;
537 assert(x.size()==y.size());
538 assert(x.size()==w.size());
546 for (i=0; i<w.size(); ++
i) {
547 if (w[i] == 0)
continue;
552 xAve /= double(nHit);
553 yAve /= double(nHit);
557 for (i=0; i<w.size(); ++
i) {
558 if (w[i] == 0)
continue;
559 Sxx += (x[
i] - xAve)*(x[i] - xAve);
560 Syy += (y[
i] - yAve)*(y[i] - yAve);
561 Sxy += (x[
i] - xAve)*(y[i] - yAve);
563 Sxx /= double(nHit-1);
564 Syy /= double(nHit-1);
565 Sxy /= double(nHit-1);
570 double xydiff = Syy - delta*Sxx;
571 double radic =
sqrt( (xydiff*xydiff) + 4*delta*Sxy*Sxy );
572 double numer = xydiff + radic;
573 double denom = 2.*Sxy;
574 M1 = ( numer ) / ( denom );
577 y1 = M1*(x1-xAve) + yAve;
578 y2 = M1*(x2-xAve) + yAve;
581 double sig2y1(0), sig2y2(0);
582 for (i=0; i<w.size(); ++
i) {
583 if (w[i] == 0)
continue;
586 dnumdx += delta*(2/double(nHit-1))*(x[i]-xAve);
587 dnumdx += 1./(2*radic)*( 2*xydiff*( -delta*2/
double(nHit-1)*(x[i]-xAve) )
588 + 8*delta*Sxy*1./double(nHit-1)*(y[
i]-yAve) );
589 dnumdy += (2/double(nHit-1))*(y[i]-yAve);
590 dnumdx += 1./(2*radic)*( 2*xydiff*( 2/
double(nHit-1)*(y[
i]-yAve) )
591 + 8*Sxy*1./double(nHit-1)*(x[
i]-xAve) );
592 double dM1dx = denom != 0. ? ( dnumdx - M1*2/double(nHit-1)*(y[
i]-yAve) ) / ( denom ) : 0.;
593 double dM1dy = denom != 0. ? ( dnumdy - M1*2/double(nHit-1)*(x[
i]-xAve) ) / ( denom ) : 0.;
604 double dy1dx = dM1dx*(x1-xAve) + M1;
605 double dy1dy = dM1dy*(x1-xAve) + 1./
double(nHit);
606 double dy2dx = dM1dx*(x2-xAve) + M1;
607 double dy2dy = dM1dy*(x2-xAve) + 1./
double(nHit);
608 sig2y1 += dy1dx*dy1dx*sigX*sigX + dy1dy*dy1dy*sigY*sigY;
609 sig2y2 += dy2dx*dy2dx*sigX*sigX + dy2dy*dy2dy*sigY*sigY;
611 sigy1 =
sqrt(sig2y1);
612 sigy2 =
sqrt(sig2y2);
635 double localHitPos[3] = { 0., 0., XY };
636 double worldHitPos[3];
638 return TVector3(worldHitPos);
646 double localLo[3] = { 0., 0., -theCell->
HalfL() };
647 double localHi[3] = { 0., 0., theCell->
HalfL() };
std::vector< double > subtrk_ErrStopY
std::string fTrackInLabel
An algorithm to perform cosmic ray track fitting.
Simple analyzer module to print to text file information required for c++ based block alignment study...
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
std::vector< double > subtrk_StartZ
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Float_t y1[n_points_granero]
void LocalToWorld(const double *local, double *world) const
unsigned short Plane() const
std::vector< double > subtrk_ErrStartZ
Float_t x1[n_points_granero]
const CellGeo * Cell(int icell) const
double DemingRegFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &sigXY, double &sigZ, double &x1, double &y1, double &x2, double &y2, double &sigy1, double &sigy2)
Find the best-fit line to a collection of points in 2-D using Deming regression.
std::vector< double > blockCenterX
Vertical planes which measure X.
std::vector< double > subtrk_ErrStopZ
A collection of associated CellHits.
const daqchannelmap::DAQChannelMap * Map() const
unsigned int MaxCell(geo::View_t view) const
A rb::Prong with full reconstructed trajectory.
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Track finder for cosmic rays.
lchan encodeLChan(int detId, plane_t plane, cell_t cell) const
An algorithm to perform cosmic ray track fitting.
Horizontal planes which measure Y.
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
unsigned short Cell() const
SplitTracks(fhicl::ParameterSet const &pset)
TH2D * xzViewSubtrkResVsX
block
print "ROW IS " print row
std::vector< double > subtrk_ErrStartX
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
std::vector< int > subtrk_Block
std::vector< double > subtrk_StopX
TH2D * yzViewSubtrkResVsZ
Collect Geo headers and supply basic geometry functions.
EDAnalyzer(Table< Config > const &config)
TH2D * yzViewSubtrkResVsY
An algorithm to perform cosmic ray track fitting.
TH2D * xzViewSubtrkResVsZ
std::vector< double > subtrk_StopY
unsigned int NYCell() const
Number of cells in the y-view.
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
std::vector< double > subtrk_ErrStartY
T * make(ARGS...args) const
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
void analyze(art::Event const &evt) override
std::vector< double > blockCenterZ
unsigned int NXCell() const
Number of cells in the x-view.
unsigned int MinCell(geo::View_t view) const
std::vector< double > subtrk_StartY
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
std::vector< double > subtrk_StopZ
assert(nhit_max >=nhit_nbins)
std::vector< double > subtrk_ErrStopX
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
virtual void InterpolateXY(double z, double &x, double &y) const
std::vector< double > blockCenterY
virtual std::vector< rb::Track > MakeTrack(const rb::Cluster *slice)=0
TVector3 Stop() const
Position of the final trajectory point.
virtual block_t computeBlock(lchan logicalchan) const =0
Which block is this lchan in?
unsigned int fMinPlanesInSubCluster
Encapsulate the cell geometry.
Track finder for cosmic rays.
std::vector< double > subtrk_StartX
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
void GetCellEndpoints(art::Ptr< rb::CellHit > cell, double *cellLo, double *cellHi)
TVector3 GeoEstimatePos(art::Ptr< rb::CellHit > cellHit, double W)