33 #include "DAQDataFormats/RawEvent.h" 34 #include "DAQDataFormats/RawTrigger.h" 35 #include "DAQDataFormats/RawTriggerMask.h" 36 #include "DAQDataFormats/RawTriggerHeader.h" 37 #include "DAQDataFormats/RawDataBlock.h" 38 #include "DAQDataFormats/RawMicroBlock.h" 39 #include "DAQDataFormats/RawMicroSlice.h" 40 #include "DAQDataFormats/RawNanoSlice.h" 42 #include "NovaDAQUtilities/TimeUtils.h" 43 #include <NovaTimingUtilities/TimingUtilities.h> 52 #include "Utilities/func/MathUtil.h" 66 #include <TTimeStamp.h> 69 #include "TTimeStamp.h" 75 #pragma link C++ class vector<double>+; 86 bool sortByNCell(Track2D &track_i, Track2D &track_j);
98 Track_Match(std::vector<Track2D>::const_iterator
x, std::vector<Track2D>::const_iterator
y):
127 void endJob()
override;
151 unsigned int totalntracksX = 0;
152 unsigned int totalntracksY = 0;
153 unsigned int total3DTracks = 0;
154 unsigned int failed3Devents = 0;
277 int nSinglemuEvents = 0;
278 int nMultimuEvents = 0;
292 double totalTracks = 0;
307 hough_peak_cut_ (p.
get<unsigned
int>(
"hough_peak_cut")),
308 track2D_min_hits_ (p.
get<unsigned
int>(
"track2D_min_hits")),
309 max_distance_between_hits_ (p.
get<double>(
"max_distance_between_hits")),
310 fake_tracks_threshold_ (p.
get<double>(
"fake_tracks_threshold")),
311 track3D_min_hits_ (p.
get<unsigned>(
"track3D_min_hits")),
312 hit_distance_cut_ (p.
get<
int>(
"hit_distance_cut")),
313 vertex_max_distance_ (p.
get<double>(
"vertex_max_distance")),
558 const int nflshits = flshits.size();
561 if (nflshits == 0) {
continue; }
564 if (TMath::Abs(particle->
PdgCode()) != 13) {
continue; }
567 const TLorentzVector& fourPosition = particle->
Position();
568 const TLorentzVector& fourMomentum = particle->
Momentum();
577 for(
unsigned int i_plane = 0, n_plane = hitsByPlaneCell.size(); i_plane < n_plane; ++i_plane) {
578 hitsByPlaneCell[i_plane].resize(numCells, 0.);
590 double entryPos[3] = {-9999, -9999, -9999};
591 double exitPos[3] = {9999, 9999, 9999};
593 double localStartPos[3] = {-9999, -9999, -9999};
594 double localEndPos[3] = {9999, 9999, 9999};
596 double entryPosTemp[3];
597 double exitPosTemp[3];
607 int FLShitsCounter = 0;
610 for (
const auto&
hit : flshits) {
612 startPlane =
hit.GetPlaneID();
613 startCell =
hit.GetCellID();
614 localStartPos[0] =
hit.GetEntryX();
615 localStartPos[1] =
hit.GetEntryY();
616 localStartPos[2] =
hit.GetEntryZ();
618 endPlane =
hit.GetPlaneID();
619 endCell =
hit.GetCellID();
620 localEndPos[0] =
hit.GetExitX();
621 localEndPos[1] =
hit.GetExitY();
622 localEndPos[2] =
hit.GetExitZ();
626 hitsByPlaneCell[
hit.GetPlaneID()][
hit.GetCellID()] += 1.;
635 if (entryPos[1] < entryPosTemp[1] || entryPos[1] < exitPosTemp[1]) {
637 entryPos[0] = entryPosTemp[0];
638 entryPos[1] = entryPosTemp[1];
639 entryPos[2] = entryPosTemp[2];
644 if (exitPos[1] > exitPosTemp[1] || exitPos[1] > entryPosTemp[1] ) {
646 exitPos[0] = entryPosTemp[0];
647 exitPos[1] = entryPosTemp[1];
648 exitPos[2] = entryPosTemp[2];
652 if (
hit.GetPlaneID() % 2) { planesX++; }
660 double pRx = exitPos[0] - entryPos[0];
661 double pRy = exitPos[1] - entryPos[1];
662 double pRz = exitPos[2] - entryPos[2];
663 double abspR =
TMath::Sqrt( pRx*pRx + pRy*pRy + pRz*pRz );
666 double deltaXsquare = (exitPos[0] - entryPos[0])*(exitPos[0] - entryPos[0]);
667 double deltaYsquare = (exitPos[1] - entryPos[1])*(exitPos[1] - entryPos[1]);
668 double deltaZsquare = (exitPos[2] - entryPos[2])*(exitPos[2] - entryPos[2]);
669 double particleLength =
TMath::Sqrt( deltaXsquare + deltaYsquare + deltaZsquare );
672 double pTheta = TMath::ACos( pRy / abspR );
673 double pZenith = TMath::Pi() - pTheta;
674 double pPhi = TMath::ATan(pRz / pRx);
677 if (pRx >= 0 && pRz >= 0) { pAzimuth = pPhi; }
678 if (pRx < 0 && pRz >= 0) { pAzimuth = TMath::Pi() + pPhi; }
679 if (pRx < 0 && pRz < 0) { pAzimuth = TMath::Pi() + pPhi; }
680 if (pRx >= 0 && pRz < 0) { pAzimuth = 2*TMath::Pi() + pPhi; }
726 particlePhi.push_back( pAzimuth * TMath::RadToDeg() );
785 unsigned int eventToCheck = 0;
791 std::cout <<
"No slices in event " << e.
event() <<
", so doing nothing\n";
799 const time_t timeSec = ts.
timeHigh();
800 struct tm* timeStruct = localtime(&timeSec);
803 unsigned int UTCyear_, UTCmonth_, UTCday_, UTChour_, UTCminute_, UTCsecond_;
805 unsigned long long int tsval = e.
time().
value();
806 const unsigned long int mask32 = 0xFFFFFFFFUL;
807 unsigned long int lup = ( tsval >> 32 ) & mask32;
808 unsigned long int llo = tsval & mask32;
809 TTimeStamp UTCts(lup, (
int)llo);
811 UTCts.GetDate(kTRUE,0,&UTCyear_, &UTCmonth_, &UTCday_);
812 UTCts.GetTime(kTRUE,0,&UTChour_, &UTCminute_, &UTCsecond_);
925 year = timeStruct->tm_year + 1900;
926 month = timeStruct->tm_mon + 1;
927 day = timeStruct->tm_mday;
928 doy = timeStruct->tm_yday + 1;
929 hour = timeStruct->tm_hour + 1;
930 minute = timeStruct->tm_min;
931 second = timeStruct->tm_sec;
937 UTCdoy = UTCts.GetDayOfYear();
1008 if (
event == eventToCheck){
1018 std::map<geo::View_t, std::vector<Track2D> > hough_tracks;
1021 for (
size_t n_slice = 0; n_slice != my_slices->size(); ++n_slice)
1023 std::vector<art::Ptr<rb::HoughResult> > hough_results;
1024 find_hough_results.
get(n_slice, hough_results);
1028 for (
auto const& hough_result : hough_results)
1031 unsigned peakInSlice = 0;
1034 for (
auto const& peak : hough_result->fPeak)
1042 double slope = -9999;
1043 double intercept = -9999;
1044 hough_result->SlopeIcept(peakInSlice++, &slope, &intercept);
1053 Track2D hough_track(view, slope, intercept);
1055 if (
event == eventToCheck) {
1057 std::cout <<
"2D: view, slope, theta: " << view <<
", " << hough_track.
slope <<
", " << TMath::Abs(TMath::ATan(hough_track.
slope)) <<
std::endl;
1067 for (
auto const&
hit : hits_planeSorted)
1070 if (
hit->View() !=
view)
continue;
1080 TVector3 hit_vector(hit_xyz);
1085 hit_Z = hit_vector.z();
1087 double wX = hit_Z - ( (hit_XY - intercept) / slope );
1088 double doca = TMath::Abs( wX * TMath::Sin(TMath::ATan(slope)) );
1091 double sinTheta = TMath::Sin(houghTheta);
1096 if (doca <= max_hit_distance) {
1117 hough_tracks[
view].push_back(hough_track);
1257 std::vector<Track2D> x_clusteredTracks;
1258 std::vector<Track2D> y_clusteredTracks;
1261 for (
auto x_track = hough_tracks[
geo::kX].begin(); x_track != hough_tracks[
geo::kX].end(); x_track++) {
1267 int trackSize = x_track->cluster.NCell();
1268 int halfTrackSize = trackSize/2;
1270 double x_slope = x_track->slope;
1271 double x_intercept = x_track->intercept;
1273 double x_theta = TMath::Abs(TMath::ATan(x_track->slope));
1274 double x_sinTheta = TMath::Sin(x_theta);
1278 if (
event == eventToCheck) {
1279 std::cout <<
"XZ: NCell, theta, halftrack, maxGap: " << trackSize <<
", " << x_theta * TMath::RadToDeg() <<
", " << halfTrackSize <<
", " << maxGap <<
std::endl;
1285 for (
int hit_i = halfTrackSize; hit_i > 1; hit_i--) {
1290 double currentHit_xyz[3];
1292 TVector3 currentHit_vector(currentHit_xyz);
1296 double currentHit_X = currentHit_vector.x();
1297 double currentHit_Z = currentHit_vector.z();
1299 double previousHit_xyz[3];
1301 TVector3 previousHit_vector(previousHit_xyz);
1303 double previousHit_X = previousHit_vector.x();
1304 double previousHit_Z = previousHit_vector.z();
1306 double deltaXsquare = (previousHit_X - currentHit_X)*(previousHit_X - currentHit_X);
1307 double deltaYsquare = (previousHit_Z - currentHit_Z)*(previousHit_Z - currentHit_Z);
1309 double distance_between_hits =
sqrt(deltaXsquare + deltaYsquare);
1314 if (distance_between_hits > maxGap) {
break; }
1320 for (
int hit_i = halfTrackSize; hit_i < trackSize - 1; hit_i++) {
1325 double currentHit_xyz[3];
1327 TVector3 currentHit_vector(currentHit_xyz);
1329 double currentHit_X = currentHit_vector.x();
1330 double currentHit_Z = currentHit_vector.z();
1332 double nextHit_xyz[3];
1334 TVector3 nextHit_vector(nextHit_xyz);
1336 double nextHit_X = nextHit_vector.x();
1337 double nextHit_Z = nextHit_vector.z();
1339 double deltaXsquare = (nextHit_X - currentHit_X)*(nextHit_X - currentHit_X);
1340 double deltaZsquare = (nextHit_Z - currentHit_Z)*(nextHit_Z - currentHit_Z);
1342 double distance_between_hits =
sqrt(deltaXsquare + deltaZsquare);
1347 if (distance_between_hits > maxGap) {
break; }
1358 for (
auto y_track = hough_tracks[
geo::kY].begin(); y_track != hough_tracks[
geo::kY].end(); y_track++) {
1364 int trackSize = y_track->cluster.NCell();
1365 int halfTrackSize = trackSize/2;
1367 double y_slope = y_track->slope;
1368 double y_intercept = y_track->intercept;
1370 double y_theta = TMath::Abs(TMath::ATan(y_track->slope));
1371 double y_sinTheta = TMath::Sin(y_theta);
1374 if (
event == eventToCheck) {
1375 std::cout <<
"YZ: NCell, theta, halftrack, maxGap: " << trackSize <<
", " << y_theta * TMath::RadToDeg() <<
", " << halfTrackSize <<
", " << maxGap <<
std::endl;
1381 for (
int hit_i = halfTrackSize; hit_i > 1; hit_i--) {
1386 double currentHit_xyz[3];
1388 TVector3 currentHit_vector(currentHit_xyz);
1390 double currentHit_Y = currentHit_vector.y();
1391 double currentHit_Z = currentHit_vector.z();
1393 double previousHit_xyz[3];
1395 TVector3 previousHit_vector(previousHit_xyz);
1397 double previousHit_Y = previousHit_vector.y();
1398 double previousHit_Z = previousHit_vector.z();
1400 double deltaYsquare = (previousHit_Y - currentHit_Y)*(previousHit_Y - currentHit_Y);
1401 double deltaZsquare = (previousHit_Z - currentHit_Z)*(previousHit_Z - currentHit_Z);
1403 double distance_between_hits =
sqrt(deltaYsquare + deltaZsquare);
1408 if (distance_between_hits > maxGap) {
break; }
1414 for (
int hit_i = halfTrackSize; hit_i < trackSize - 1; hit_i++) {
1419 double currentHit_xyz[3];
1421 TVector3 currentHit_vector(currentHit_xyz);
1423 double currentHit_Y = currentHit_vector.y();
1424 double currentHit_Z = currentHit_vector.z();
1426 double nextHit_xyz[3];
1428 TVector3 nextHit_vector(nextHit_xyz);
1430 double nextHit_Y = nextHit_vector.y();
1431 double nextHit_Z = nextHit_vector.z();
1433 double deltaYsquare = (nextHit_Y - currentHit_Y)*(nextHit_Y - currentHit_Y);
1434 double deltaZsquare = (nextHit_Z - currentHit_Z)*(nextHit_Z - currentHit_Z);
1436 double distance_between_hits =
sqrt(deltaYsquare + deltaZsquare);
1441 if (distance_between_hits > maxGap) {
break; }
1457 if (x_clusteredTracks.size() == 0 || y_clusteredTracks.size() == 0) {
1565 std::sort(x_clusteredTracks.begin(), x_clusteredTracks.end(),
htk::sortByNCell);
1566 std::sort(y_clusteredTracks.begin(), y_clusteredTracks.end(),
htk::sortByNCell);
1569 std::reverse(x_clusteredTracks.begin(), x_clusteredTracks.end());
1570 std::reverse(y_clusteredTracks.begin(), y_clusteredTracks.end());
1574 std::vector<Track2D> x_goodTracks;
1575 std::vector<Track2D> y_goodTracks;
1577 if (
event == eventToCheck){
1580 std::cerr <<
"x_clusteredTrack, y_clusteredTrack sizes: " << x_clusteredTracks.size() <<
", " << y_clusteredTracks.size() <<
std::endl;
1591 bool **isGood2DxTrackMatrix;
1592 isGood2DxTrackMatrix =
new bool*[x_clusteredTracks.size()];
1594 for (
unsigned int p = 0;
p < x_clusteredTracks.size();
p++) {
1596 isGood2DxTrackMatrix[
p] =
new bool[x_clusteredTracks.size()];
1601 for (
unsigned int n = 0;
n < x_clusteredTracks.size();
n++) {
1603 for (
unsigned int m = 0;
m < x_clusteredTracks.size();
m++) {
1605 isGood2DxTrackMatrix[
n][
m] =
true;
1611 for (
unsigned int i = 0;
i < x_clusteredTracks.size() - 1;
i++)
1614 Track2D x_track_i = x_clusteredTracks[
i];
1617 for (
unsigned int j =
i + 1;
j < x_clusteredTracks.size();
j++)
1620 Track2D x_track_j = x_clusteredTracks[
j];
1623 for (
auto const &hit_i: x_hits_i)
1626 double hit_i_xyz[3];
1628 TVector3 hit_i_vector(hit_i_xyz);
1630 double hit_i_Z = hit_i_vector.z();
1632 for (
auto const &hit_j: x_hits_j)
1635 double hit_j_xyz[3];
1637 TVector3 hit_j_vector(hit_j_xyz);
1639 double hit_j_Z = hit_j_vector.z();
1643 isGood2DxTrackMatrix[
i][
j] =
false;
1649 if (isGood2DxTrackMatrix[
i][
j] ==
false) {
continue; }
1656 for (
unsigned int n = 0;
n < x_clusteredTracks.size();
n++) {
1658 isGood2DTrack =
true;
1660 for (
unsigned int m = 0;
m < x_clusteredTracks.size();
m++) {
1662 if (isGood2DxTrackMatrix[
m][
n] ==
false) { isGood2DTrack =
false; }
1666 if (isGood2DTrack ==
true) { x_goodTracks.push_back(x_clusteredTracks[
n]); }
1675 bool **isGood2DyTrackMatrix;
1676 isGood2DyTrackMatrix =
new bool*[y_clusteredTracks.size()];
1678 for (
unsigned int p = 0;
p < y_clusteredTracks.size();
p++) {
1680 isGood2DyTrackMatrix[
p] =
new bool[y_clusteredTracks.size()];
1683 for (
unsigned int n = 0;
n < y_clusteredTracks.size();
n++) {
1685 for (
unsigned int m = 0;
m < y_clusteredTracks.size();
m++) {
1687 isGood2DyTrackMatrix[
n][
m] =
true;
1693 for (
unsigned int i = 0;
i < y_clusteredTracks.size() - 1;
i++)
1696 Track2D y_track_i = y_clusteredTracks[
i];
1699 for (
unsigned int j =
i + 1;
j < y_clusteredTracks.size();
j++)
1702 Track2D y_track_j = y_clusteredTracks[
j];
1705 for (
auto const &hit_i: y_hits_i)
1708 double hit_i_xyz[3];
1710 TVector3 hit_i_vector(hit_i_xyz);
1712 double hit_i_Z = hit_i_vector.z();
1714 for (
auto const &hit_j: y_hits_j)
1717 double hit_j_xyz[3];
1719 TVector3 hit_j_vector(hit_j_xyz);
1721 double hit_j_Z = hit_j_vector.z();
1725 isGood2DyTrackMatrix[
i][
j] =
false;
1731 if (isGood2DyTrackMatrix[
i][
j] ==
false) {
continue; }
1738 for (
unsigned int n = 0;
n < y_clusteredTracks.size();
n++) {
1740 isGood2DTrack =
true;
1742 for (
unsigned int m = 0;
m < y_clusteredTracks.size();
m++) {
1744 if (isGood2DyTrackMatrix[
m][
n] ==
false) { isGood2DTrack =
false; }
1748 if (isGood2DTrack ==
true) { y_goodTracks.push_back(y_clusteredTracks[
n]); }
1757 if (
event == eventToCheck) {
1765 for (
unsigned int n = 0;
n < x_clusteredTracks.size();
n++) {
1773 for (
unsigned int n = 0;
n < x_clusteredTracks.size();
n++) {
1777 for (
unsigned int m = 0;
m < x_clusteredTracks.size();
m++) {
1779 if (
m == 0) {
std::cerr << isGood2DxTrackMatrix[
n][
m]; }
1780 if (
m > 0) {
std::cerr <<
" " << isGood2DxTrackMatrix[
n][
m]; }
1797 for (
unsigned int n = 0;
n < y_clusteredTracks.size();
n++) {
1805 for (
unsigned int n = 0;
n < y_clusteredTracks.size();
n++) {
1809 for (
unsigned int m = 0;
m < y_clusteredTracks.size();
m++) {
1811 if (
m == 0) {
std::cerr << isGood2DyTrackMatrix[
n][
m]; }
1812 if (
m > 0) {
std::cerr <<
" " << isGood2DyTrackMatrix[
n][
m]; }
1828 double x_hit_startX = 9999;
1829 double x_hit_startZ = 9999;
1832 int xz_startPlane = 0;
1833 int xz_endPlane = 0;
1834 double xz_hit_start_time = 9999;
1837 double y_hit_startY = 9999;
1838 double y_hit_startZ = 9999;
1841 int yz_startPlane = 0;
1842 int yz_endPlane = 0;
1843 double yz_hit_start_time = 9999;
1846 double x_hit_endX = -9999;
1847 double x_hit_endZ = -9999;
1849 double xz_hit_end_time = 9999;
1851 double y_hit_endY = -9999;
1852 double y_hit_endZ = -9999;
1854 double yz_hit_end_time = 9999;
1861 std::reverse(x_goodTracks.begin(), x_goodTracks.end());
1862 std::reverse(y_goodTracks.begin(), y_goodTracks.end());
1866 std::vector<Track_Match> track_matches;
1867 std::unique_ptr<std::vector<Track3D> > tracks_3D(
new std::vector<Track3D>);
1872 for (
auto x_track = x_goodTracks.begin(); x_track != x_goodTracks.end(); x_track++)
1879 double x_track_meanTime = 0;
1880 unsigned int nhitsXZ = 0;
1881 int this_xz_plane = -99;
1884 for (
auto const& x_hit : x_hits_planeSorted)
1888 double x_hit_xyz[3];
1890 TVector3 x_hit_vector(x_hit_xyz);
1894 x_hit_startX = x_hit_vector.x();
1895 x_hit_startZ = x_hit_vector.z();
1897 xz_startPlane = x_hit->Plane();
1898 xz_hit_start_time = x_hit->TNS();
1902 if (nhitsXZ == x_track->cluster.NCell() - 1) {
1904 x_hit_endX = x_hit_vector.x();
1905 x_hit_endZ = x_hit_vector.z();
1907 xz_endPlane = x_hit->Plane();
1908 xz_hit_end_time = x_hit->TNS();
1912 if (x_hit->Plane() != this_xz_plane) {
1914 this_xz_plane = x_hit->Plane();
1919 x_track_meanTime = x_track_meanTime + x_hit->TNS();
1927 for (
auto y_track = y_goodTracks.begin(); y_track != y_goodTracks.end(); y_track++)
1933 double y_track_meanTime = 0;
1934 unsigned int nhitsYZ = 0;
1935 int this_yz_plane = -99;
1938 for (
auto const& y_hit : y_hits_planeSorted)
1942 double y_hit_xyz[3];
1944 TVector3 y_hit_vector(y_hit_xyz);
1948 y_hit_startY = y_hit_vector.y();
1949 y_hit_startZ = y_hit_vector.z();
1951 yz_startPlane = y_hit->Plane();
1952 yz_hit_start_time = y_hit->TNS();
1957 if (nhitsYZ == y_track->cluster.NCell() - 1) {
1959 y_hit_endY = y_hit_vector.y();
1960 y_hit_endZ = y_hit_vector.z();
1962 yz_endPlane = y_hit->Plane();
1963 yz_hit_end_time = y_hit->TNS();
1968 if (y_hit->Plane() != this_yz_plane) {
1970 this_yz_plane = y_hit->Plane();
1975 y_track_meanTime = y_track_meanTime + y_hit->TNS();
1980 x_track_meanTime = x_track_meanTime / x_track->cluster.NCell();
1981 y_track_meanTime = y_track_meanTime / y_track->cluster.NCell();
1989 if (
event == eventToCheck){
1993 std::cerr <<
"NCell : " << x_track->cluster.NCell() <<
", " << y_track->cluster.NCell() <<
std::endl;
1994 std::cerr <<
"Start : (" << x_hit_startX <<
", " << x_hit_startZ <<
"), (" << y_hit_startY <<
", " << y_hit_startZ <<
")" <<
std::endl;
1995 std::cerr <<
"End : (" << x_hit_endX <<
", " << x_hit_endZ <<
"), (" << y_hit_endY <<
", " << y_hit_endZ <<
")" <<
std::endl;
1996 std::cerr <<
"Theta : " << TMath::ATan(x_track->slope) * TMath::RadToDeg() <<
", " << TMath::ATan(y_track->slope) * TMath::RadToDeg() <<
std::endl;
2007 bool isGood3DTrack =
true;
2011 for (
unsigned int i = 0;
i <
startZ.size();
i++) {
2014 if (
thetaXZ.at(
i) == (TMath::ATan(x_track->slope) * TMath::RadToDeg()) ||
thetaYZ.at(
i) == (TMath::ATan(y_track->slope) * TMath::RadToDeg()) ) {
2016 isGood3DTrack =
false;
2023 if (isGood3DTrack ==
true) {
2028 if (y_hit_startY < y_hit_endY) {
2033 hit_temp = y_hit_startY;
2034 y_hit_startY = y_hit_endY;
2035 y_hit_endY = hit_temp;
2038 hit_temp = y_hit_startZ;
2039 y_hit_startZ = y_hit_endZ;
2040 y_hit_endZ = hit_temp;
2043 hit_temp = yz_hit_start_time;
2044 yz_hit_start_time = yz_hit_end_time;
2045 yz_hit_end_time = hit_temp;
2048 if ( (y_hit_startZ - y_hit_endZ) > 0 && (x_hit_startZ - x_hit_endZ) < 0 ) {
2051 hit_temp = x_hit_startX;
2052 x_hit_startX = x_hit_endX;
2053 x_hit_endX = hit_temp;
2056 hit_temp = x_hit_startZ;
2057 x_hit_startZ = x_hit_endZ;
2058 x_hit_endZ = hit_temp;
2061 hit_temp = xz_hit_start_time;
2062 xz_hit_start_time = xz_hit_end_time;
2063 xz_hit_end_time = hit_temp;
2065 hit_temp = xz_startPlane;
2066 xz_startPlane = xz_endPlane;
2067 xz_endPlane = hit_temp;
2071 if ( (y_hit_startZ - y_hit_endZ) < 0 && (x_hit_startZ - x_hit_endZ) > 0 ) {
2074 hit_temp = x_hit_startX;
2075 x_hit_startX = x_hit_endX;
2076 x_hit_endX = hit_temp;
2079 hit_temp = x_hit_startZ;
2080 x_hit_startZ = x_hit_endZ;
2081 x_hit_endZ = hit_temp;
2084 hit_temp = xz_hit_start_time;
2085 xz_hit_start_time = xz_hit_end_time;
2086 xz_hit_end_time = hit_temp;
2088 hit_temp = xz_startPlane;
2089 xz_startPlane = xz_endPlane;
2090 xz_endPlane = hit_temp;
2098 double Rx = x_hit_endX - x_hit_startX;
2099 double Ry = y_hit_endY - y_hit_startY;
2100 double Rz = y_hit_endZ - y_hit_startZ;
2101 double absR =
TMath::Sqrt( Rx*Rx + Ry*Ry + Rz*Rz );
2103 double trackTheta = TMath::ACos( Ry / absR );
2104 double zen = TMath::Pi() - trackTheta;
2105 double trackPhi = TMath::ATan(Rz / Rx);
2108 if (Rx > 0 && Rz > 0) { azi = trackPhi; }
2109 if (Rx < 0 && Rz > 0) { azi = TMath::Pi() + trackPhi; }
2110 if (Rx < 0 && Rz < 0) { azi = TMath::Pi() + trackPhi; }
2111 if (Rx > 0 && Rz < 0) { azi = 2*TMath::Pi() + trackPhi; }
2112 if (Rx == 0 && Rz > 0) { azi = TMath::Pi()/2; }
2113 if (Rx == 0 && Rz < 0) { azi = 3*TMath::Pi()/2; }
2114 if (Rx > 0 && Rz == 0) { azi = 0; }
2115 if (Rx < 0 && Rz == 0) { azi = TMath::Pi(); }
2117 Track3D track3D(*x_track, *y_track);
2119 double track3DmeanTime = (x_track_meanTime + y_track_meanTime) / 2;
2125 nplanes.push_back(xz_nplanes + yz_nplanes);
2131 len.push_back(absR);
2136 if (xz_hit_start_time < yz_hit_start_time) {
startTime.push_back(xz_hit_start_time); }
2137 else {
startTime.push_back(yz_hit_start_time); }
2138 meanTime.push_back(track3DmeanTime);
2142 startX.push_back(x_hit_startX);
2143 startY.push_back(y_hit_startY);
2144 if (x_hit_startZ > y_hit_startZ)
startZ.push_back(x_hit_startZ);
2145 else startZ.push_back(y_hit_startZ);
2148 endX.push_back(x_hit_endX);
2149 endY.push_back(y_hit_endY);
2150 if (x_hit_endZ > y_hit_endZ)
endZ.push_back(x_hit_endZ);
2151 else endZ.push_back(y_hit_endZ);
2153 dirX.push_back( Rx/absR );
2154 dirY.push_back( Ry/absR );
2155 dirZ.push_back( Rz/absR );
2157 thetaXZ.push_back( TMath::ATan(x_track->slope) * TMath::RadToDeg() );
2158 thetaYZ.push_back( TMath::ATan(y_track->slope) * TMath::RadToDeg() );
2160 theta.push_back( zen * TMath::RadToDeg() );
2161 phi.push_back( azi * TMath::RadToDeg() );
2162 cosTheta.push_back( TMath::Cos(zen) );
2163 cosPhi.push_back( TMath::Cos(azi) );
2166 double trueNorthCorrection = 0.668112;
2168 zenith.push_back( zen * TMath::RadToDeg() );
2169 if (azi - trueNorthCorrection > 0) {
azimuth.push_back( (azi - trueNorthCorrection) * TMath::RadToDeg() ); }
2170 if (azi - trueNorthCorrection < 0) {
azimuth.push_back( (TMath::Pi()*2 - TMath::Abs(azi - trueNorthCorrection)) * TMath::RadToDeg() ); }
2172 cosAzimuth.push_back( TMath::Cos(azi - trueNorthCorrection));
2175 track_matches.emplace_back(x_track, y_track);
2176 tracks_3D->push_back(track3D);
2178 if (
event == eventToCheck){
2183 std::cerr <<
"NCell : " << x_track->cluster.NCell() <<
", " << y_track->cluster.NCell() <<
std::endl;
2184 std::cerr <<
"Start : (" << x_hit_startX <<
", " << x_hit_startZ <<
"), (" << y_hit_startY <<
", " << y_hit_startZ <<
")" <<
std::endl;
2185 std::cerr <<
"End : (" << x_hit_endX <<
", " << x_hit_endZ <<
"), (" << y_hit_endY <<
", " << y_hit_endZ <<
")" <<
std::endl;
2186 std::cerr <<
"Theta : " << TMath::ATan(x_track->slope) * TMath::RadToDeg() <<
", " << TMath::ATan(y_track->slope) * TMath::RadToDeg() <<
std::endl;
2300 for (
unsigned int k = 0; k < x_clusteredTracks.size(); k++) {
2302 delete[] isGood2DxTrackMatrix[k];
2306 delete[] isGood2DxTrackMatrix;
2307 isGood2DxTrackMatrix = 0;
2309 for (
unsigned int k = 0; k < y_clusteredTracks.size(); k++) {
2311 delete[] isGood2DyTrackMatrix[k];
2315 delete[] isGood2DyTrackMatrix;
2316 isGood2DyTrackMatrix = 0;
unsigned long long triggerStartUnixTime
std::vector< double > particleTheta
SubRunNumber_t subRun() const
std::vector< double > nplanesYZ
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
std::vector< double > calE
std::vector< double > particleEndX
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
std::vector< double > cosAzimuth
std::vector< double > particleDirX
std::vector< double > particleNplanesXZ
std::vector< double > particleStartY
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
std::vector< double > endZ
std::vector< Track2D >::const_iterator x_track
std::vector< double > cosZenith
std::vector< double > startZ
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
void LocalToWorld(const double *local, double *world) const
unsigned short Plane() const
std::vector< double > endPlaneYZ
std::vector< double > dirZ
rb::Cluster cluster() const
const CellGeo * Cell(int icell) const
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
constexpr std::uint32_t timeHigh() const
list_type::const_iterator const_iterator
std::vector< double > particleCosPhi
Vertical planes which measure X.
std::vector< double > endY
unsigned int Ncells() const
Number of cells in this plane.
std::vector< double > particlePx
std::vector< double > cosTheta
const PlaneGeo * Plane(unsigned int i) const
std::vector< double > particleEndZ
std::vector< double > particleE
DEFINE_ART_MODULE(TestTMapFile)
std::vector< double > meanTime
unsigned track2D_min_hits_
bool convertNovaTimeToUnixTime(uint64_t const &inputNovaTime, struct timespec &outputUnixTime)
std::vector< double > endX
constexpr TimeValue_t value() const
Horizontal planes which measure Y.
std::vector< double > startX
bool sortByFLSHits(sim::Particle particle_i, sim::Particle particle_j)
const Var kY([](const caf::SRProxy *sr){float tmp=0.f;if(sr->mc.nu.empty()) return tmp;tmp=sr->mc.nu[0].y;return tmp;})
std::vector< double > particleNplanesYZ
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
std::vector< double > particleEndY
std::vector< double > startPlaneXZ
std::vector< double > startTime
unsigned short Cell() const
unsigned int failed3Devents
unsigned int totalntracksX
unsigned int totalntracksY
std::vector< double > particleMass
std::vector< double > particleFLSHits
std::vector< double > nplanes
art::ServiceHandle< cheat::BackTracker > bt
std::vector< double > houghIntercept
unsigned int hough_peak_cut_
std::vector< double > houghView
Collect Geo headers and supply basic geometry functions.
std::vector< double > particleDirY
std::vector< double > houghSlope
double vertex_max_distance_
bool sortByNCell(Track2D &track_i, Track2D &track_j)
Track_Match(std::vector< Track2D >::const_iterator x, std::vector< Track2D >::const_iterator y)
std::vector< double > thetaYZ
std::vector< double > particleT
EventNumber_t event() const
unsigned long long timestamp
std::vector< double > len
void analyze(art::Event const &e) override
std::vector< double > particleDirZ
std::vector< double > houghHeight
T * make(ARGS...args) const
Data resulting from a Hough transform on the cell hit positions.
std::vector< double > nplanesXZ
std::vector< double > startY
std::vector< double > particlePhi
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Var Sqrt(const Var &v)
Use to take sqrt of a var.
unsigned long long fTriggerTimingMarker_TimeStart
std::vector< double > thetaXZ
std::vector< double > particleStartZ
std::vector< double > phi
const std::set< unsigned int > & GetPlanesByView(View_t v=kXorY) const
const TLorentzVector & Momentum(const int i=0) const
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
std::vector< double > particlePy
HoughTrack(fhicl::ParameterSet const &p)
std::vector< double > theta
std::vector< double > zenith
std::vector< double > nhits
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
std::vector< double > particleLen
std::vector< double > cosPhi
unsigned int total3DTracks
size_type get(size_type i, reference item, data_reference data) const
art::ServiceHandle< geo::Geometry > geometry_
unsigned track3D_min_hits_
double max_distance_between_hits_
double fake_tracks_threshold_
std::vector< double > particleStartX
double startTimeInSeconds
std::vector< double > dirY
std::vector< double > startPlaneYZ
std::vector< double > particlePDG
std::vector< double > houghTheta
Encapsulate the geometry of one entire detector (near, far, ndos)
std::vector< double > azimuth
std::vector< double > houghRho
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
std::vector< Track2D >::const_iterator y_track
std::vector< double > particleCosTheta
std::vector< double > endPlaneXZ
std::vector< double > particlePz
std::vector< double > dirX