14 #include <boost/format.hpp> 32 #include "TTimeStamp.h" 62 #include "Utilities/AssociationUtil.h" 63 #include "Utilities/func/MathUtil.h" 64 #include "NovaDAQConventions/DAQConventions.h" 90 bool DoTracksIntersectYZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2)
const;
92 bool DoTracksIntersectXZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2)
const;
96 std::map<std::string, TH1*>
h_;
104 void z2north (std::vector<TVector3> track_i,
105 std::vector<TVector3> track_o);
134 void AngleBetweentrk (std::vector<TVector3> &trkcol, TVector3 &trkvec,
double &
sigma, std::vector<double> &sigmacol);
297 double phideg = 332.0 + (03.0/60.0) + (58.071769/3600);
308 std::vector<unsigned>
Geventcol = {6930,6956, 6983, 6990, 7006, 7022, 7083, 7090, 7089, 7103,
309 7111, 7133, 7151, 7177, 7182, 7216, 7241, 7337, 7392,
310 7485, 7500, 7536, 7555, 7567, 7714, 7728, 7746, 7748,
311 7760, 7766, 7796, 7864, 7934, 7962, 7980, 7981, 8002,
312 8050, 8072, 8093, 8094, 8156, 8158, 8235, 8277, 8308,
313 8366, 8411, 8645, 8489, 8427, 8443, 8456, 8527, 8545,
314 8560, 8564, 8633, 8643, 8695, 8699, 8727, 8748, 8953,
315 9051, 9081, 9086, 9106, 9137, 9167, 9168};
317 std::vector<unsigned>
BadEventcol = {7392, 7702, 7672, 7839, 7629, 7587, 8149, 8065};
347 std::map<std::string, double>
srt_;
350 std::map<std::string, double>
et_;
387 std::vector<std::string> subrun_names {
"exposure",
"totevents",
"tottracks",
"passedtracks",
"passedevents",
"cosmicrate"};
388 for (
auto const&
name : subrun_names)
392 std::vector<std::string> event_names {
"time",
"timestamp"};
393 for (
auto const&
name : event_names)
399 "Closest distance between two tracks; Distance (m) ;N 3DTracks / 1 m",
403 "enddist - startdist (between 2 tracks) ; Distance (cm) ;N 3DTracks / 2.0 cm",
408 "#Delta #theta ;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees",
409 9000,-0.005, 89.995);
412 "#Delta #theta ;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees",
413 9000,-0.005, 89.995);
416 "#Delta #theta;#Delta #theta (degrees);N 3DTracks / 0.01 degrees",
417 9000,-0.005, 89.995);
421 "#Delta #theta;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees",
422 9000,-0.005, 89.995);
425 "#Delta #theta;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees",
426 9000,-0.005, 89.995);
429 "Azimutal Angle;Azimuth (degrees); N Tracks / 1.0^{o}",
433 "Altitude Angle;Altitude (degrees);N Tracks / 1.0^{o}",
438 "(#Delta#theta)^{2} /2 ; (#Delta#theta)^{2} /2 (deg^{2}); N Muon Pairs/ 0.01 (deg^{2})",
442 "#Delta #theta_{average} vs Tracks per Slice; #Delta #theta_{average} per slice (Degrees); N Tracks per slice",
448 "Azimuth_{average} vs Tracks per Slice; Azimuth_{average} per slice (Degrees); N Tracks per slice",
454 "#sigma_{average} vs Altitude_{average} per slice; #sigma_{average} per slice (Degrees); Altitude_{average} per slice",
461 "Tracks Per Event;N 3D Tracks; NEvents",
465 "Tracks Per Slice;N 3D Tracks; NSlices",
469 "3D Tracks Passed Per Slice;N 3D Tracks; NSlices",
473 "N Slices with PASSED Tracks per Event; N Slices per Event; N Events",
478 fExposure = tfs->
make<TH1D>(
"exposure",
";Exposure (s)", 1000, 0., 10);
481 fCosmicRate = tfs->
make<TH1D>(
"cosmicrate",
";Cosmic Rate (kHz)", 500, 0., 100);
484 tottracks = tfs->
make<TH1D>(
"alltrkcount",
"ntracks;adc/track", 100, -0.5, 999.5);
485 conttracks = tfs->
make<TH1D>(
"conttrkcount",
"ntracks;adc/track", 100, -0.5, 999.5);
486 inttracks = tfs->
make<TH1D>(
"inttrkcount",
"ntracks;adc/track", 100, -0.5, 999.5);
487 opptracks = tfs->
make<TH1D>(
"opptrkcount",
"ntracks;adc/track", 100, -0.5, 999.5);
488 lentracks = tfs->
make<TH1D>(
"lentrkcount",
"ntracks;adc/track", 100, -0.5, 999.5);
498 if(filtered && tracks.size() > 0){
500 mf::LogVerbatim(
"KalmanTrackMerge")<<
"Filtered slice has non-zero number of reconstructed tracks. This should never happen.";
507 if(slice->
IsNoise() && tracks.size() > 0){
509 mf::LogVerbatim(
"KalmanTrackMerge")<<
"Noise slice has non-zero number of reconstructed tracks. This should never happen.";
522 for(
unsigned int itrack = 0; itrack < tracks.size(); ++itrack){
528 for(
unsigned int itrhit = 0; itrhit < trackhits.
size(); ++itrhit){
530 bool foundhit =
false;
534 for(
unsigned int islhit = 0; islhit < slicehits.
size(); ++islhit){
535 if(trackhits[itrhit] == slicehits[islhit]){
549 mf::LogVerbatim(
"KalmanTrackAna")<<
"Track made from hits on a slice other than the slice the track is associated to. This should never happen.";
571 TVector3 highest, lowest;
584 TVector3
length = highest - lowest;
602 unitvec_newxz.SetXYZ( (unitvec.X()*
cos(
phi) ) - (unitvec.Z()*
sin(
phi) ),
615 for (
unsigned int u = 0;
u < trkcol.size();
u++){
619 avgtrk.SetXYZ(sumtrk.X()/trkcol.size(),
620 sumtrk.Y()/trkcol.size(),
621 sumtrk.Z()/trkcol.size());
668 TVector3 highest, lowest;
681 TVector3
length = highest - lowest;
701 for (
unsigned int u = 0;
u < pramcol.size();
u++){
704 avgpram = sum/pramcol.size();
711 double scalarproduct;
712 scalarproduct = (
avgtrk.X()*trkvec.X()) + (
avgtrk.Y()*trkvec.Y()) + (
avgtrk.Z()*trkvec.Z());
713 sigma = (
acos(scalarproduct))*180/
M_PI;
722 double scalarproduct;
723 for (
unsigned int n = 0;
n < trkcol.size();
n++) {
725 if (trkcol[
n] == trkvec)
continue;
726 scalarproduct = (trkcol[
n].X()*trkvec.X()) + (trkcol[
n].
Y()*trkvec.Y()) + (trkcol[
n].
Z()*trkvec.Z());
727 sigma = (
acos(scalarproduct))*180/
M_PI;
728 sigmacol.push_back(sigma);
737 unsigned int sametrk = 0;
738 for (
unsigned int i = 0;
i < trkcollection.size() ; ++
i)
743 if (ftrack == trackin )
763 uxz[0] = dirvec1.X();
764 uxz[1] = dirvec1.Z();
767 vxz[0] = dirvec2.X();
768 vxz[1] = dirvec2.Z();
771 wxz[0] = highest1.X() - highest2.X();
772 wxz[1] = highest1.Z() - highest2.Z();
774 double I0xz[2] = {0.};
776 double Dxz = uxz[0]*vxz[1] - uxz[1]*vxz[0];
778 double highxz1[2] = {highest1.X(), highest1.Z()};
783 double sIxz = (vxz[0]*wxz[1] - vxz[1]*wxz[0]) / Dxz;
786 double tIxz = (uxz[0]*wxz[1] - uxz[1]*wxz[0]) / Dxz;
788 if ( (sIxz < 0 || sIxz > 1) &&
789 (tIxz < 0 || tIxz > 1) )
797 else if ( (sIxz >= 0 && sIxz <= 1) &&
798 (tIxz >= 0 && tIxz <= 1) ){
800 I0xz[0] = highxz1[0] + sIxz * uxz[0];
801 I0xz[1] = highxz1[1] + sIxz * uxz[1];
806 if (TMath::Abs(I0xz[0]) < 800) {
807 std::cout <<
"There is an XZ intersect on track AT point (" << I0xz[0] <<
", 0, " << I0xz[1] <<
")." <<
std::endl;
825 uyz[0] = dirvec1.Y();
826 uyz[1] = dirvec1.Z();
829 vyz[0] = dirvec2.Y();
830 vyz[1] = dirvec2.Z();
833 wyz[0] = highest1.Y() - highest2.Y();
834 wyz[1] = highest1.Z() - highest2.Z();
836 double highyz1[2] = {highest1.Y(), highest1.Z()};
838 double I0yz[2] = {0.};
840 double Dyz = uyz[0]*vyz[1] - uyz[1]*vyz[0];
845 double sIyz = (vyz[0]*wyz[1] - vyz[1]*wyz[0]) / Dyz;
848 double tIyz = (uyz[0]*wyz[1] - uyz[1]*wyz[0]) / Dyz;
850 if ( (sIyz < 0 || sIyz > 1) &&
851 (tIyz < 0 || tIyz > 1) )
859 else if ( (sIyz >= 0 && sIyz <= 1) &&
860 (tIyz >= 0 && tIyz <= 1) ) {
862 I0yz[0] = highyz1[0] + sIyz * uyz[0];
863 I0yz[1] = highyz1[1] + sIyz * uyz[1];
870 if (TMath::Abs(I0yz[0]) < 800 ) {
871 std::cout <<
"There is a YZ intersect on track AT point (0," << I0yz[0] <<
", " << I0yz[1] <<
")." <<
std::endl;
890 for (
unsigned int j = i+1;
j < (tracks.size() - 1); ++
j)
894 TVector3 recoStart1, recoStart2, recoEnd1, recoEnd2;
897 recoStart1.SetXYZ(ftrack->
Start().X(),
899 ftrack->
Start().Z());
901 recoEnd1.SetXYZ(ftrack->
Stop().X(),
907 TVector3 highest1, lowest1;
908 if (recoStart1.Y() > recoEnd1.Y())
910 highest1 = recoStart1;
919 recoStart2.SetXYZ(track->
Start().X(),
923 recoEnd2.SetXYZ(track->
Stop().X(),
927 TVector3 highest2, lowest2;
928 if (recoStart2.Y() > recoEnd2.Y())
930 highest2 = recoStart2;
939 TVector3 dirvec1, dirvec2, diffvec;
941 dirvec1 = highest1-lowest1;
942 dirvec2 = highest2 - lowest2;
944 double slope1Y = (highest1.Y() - lowest1.Y())/(highest1.Z() - lowest1.Z());
946 double slope2Y = (highest2.Y() - lowest2.Y())/(highest2.Z() - lowest2.Z());
949 double slope1X = (highest1.X() - lowest1.X())/(highest1.Z() - lowest1.Z());
951 double slope2X = (highest2.X() - lowest2.X())/(highest2.Z() - lowest2.Z());
955 double signtestX = slope1X*slope2X;
957 double signtestY = slope1Y*slope2Y;
960 if (signtestX < 0 || signtestY < 0)
963 std::cout <<
"Track " << i <<
" and Track " <<
j <<
" have opposite slopes. " <<
std::endl;
996 for(
unsigned int i = 0;
i<slicecol->size();++
i)
1016 TTimeStamp nowevent = tts;
1018 et_[
"time"] = tts.AsDouble();
1019 et_[
"timestamp"] = tts;
1039 outfile <<
" *************************************************************** " <<
std::endl;
1086 for (
unsigned int islice = 0; islice < slicelist.
size(); ++islice)
1089 std::vector<art::Ptr<rb::Track> > tracks;
1092 if(fmtrack.isValid())
1094 tracks = fmtrack.at(islice);
1101 std::vector<art::Ptr<rb::Track> > n3Dtracks;
1102 std::vector<art::Ptr<rb::Track> > passedtracks;
1104 std::vector<art::Ptr<rb::Track> > nocrosstracks;
1106 std::vector<art::Ptr<rb::Track> > sameslopetracks;
1108 std::vector<art::Ptr<rb::Track > > minlengthtracks;
1110 double scalarproduct;
1114 std::vector<TVector3> alltracksvec;
1115 std::vector<TVector3> unittracks;
1116 std::vector<TVector3> passedveccol;
1117 std::vector<TVector3> samesignveccol;
1118 std::vector<TVector3> nocrossveccol;
1119 std::vector<TVector3> minlengthveccol;
1120 std::vector<TVector3> sigmapassedveccol;
1121 std::vector<art::Ptr<rb::Track> > sigmapassedtracks;
1124 if(slicelist[islice]->IsNoise())
continue;
1151 for (
unsigned int i = 0;
i < tracks.size(); ++
i)
1161 z2north(unitvec, unitvec_newxz);
1165 alltracksvec.push_back(
trkvec);
1166 n3Dtracks.push_back(track);
1174 if(n3Dtracks.size() < 2)
continue;
1176 fnAll += n3Dtracks.size();
1179 for (
unsigned int i=0;
i< n3Dtracks.size(); ++
i)
1187 z2north(unitvec, unitvec_newxz);
1193 std::vector<double> sigmacol;
1195 for (
unsigned int j =
i+1;
j < (
trkcol.size() - 1);
j++) {
1197 deltatheta = (
acos(scalarproduct))*180/
M_PI;
1207 for (
unsigned int i=0;
i<n3Dtracks.size();
i++)
1220 double xyz[3] = { track->
Stop().X(), track->
Stop().Y(), track->
Stop().Z() };
1221 double xyzstart[3] = {track->
Start().X(), track->
Start().Y(), track->
Start().Z() };
1224 double highest1[3], lowest1[3];
1225 if (xyzstart[1] > xyz[1])
1227 highest1[0] = xyzstart[0];
1228 highest1[1] = xyzstart[1];
1229 highest1[2] = xyzstart[2];
1230 lowest1[0] = xyz[0];
1231 lowest1[1] = xyz[1];
1232 lowest1[2] = xyz[2];
1239 lowest1[0]=xyzstart[0];
1240 lowest1[1]=xyzstart[1];
1241 lowest1[2]=xyzstart[2];
1277 if (dist2edgeStart > 100.0){
1280 if (dist2edgeOLD > 100.0){
1287 passedtracks.push_back(track);
1297 if (passedtracks.size() < 2)
continue;
1298 fnCont += passedtracks.size();
1299 for (
unsigned int i=0;
i<passedtracks.size(); ++
i)
1308 z2north(unitvec, unitvec_newxz);
1312 passedveccol.push_back(
trkvec);
1320 for (
unsigned int i=0;
i<passedtracks.size(); ++
i)
1328 z2north(unitvec, unitvec_newxz);
1333 std::vector<double> sigmacol;
1335 for (
unsigned int j =
i+1;
j < (
trkcol.size() - 1);
j++) {
1337 deltatheta = (
acos(scalarproduct))*180/
M_PI;
1350 for (
unsigned int i = 0;
i < passedtracks.size(); ++
i)
1355 unsigned int XZcrosscount = 0;
1356 unsigned int YZcrosscount = 0;
1358 unsigned int noXZcross = 0;
1359 unsigned int noYZcross = 0;
1361 for (
unsigned int j =
i + 1;
j < (passedtracks.size() - 1); ++
j)
1365 std::vector<art::Ptr<rb::Track> > trkcollection = passedtracks;
1367 TVector3 recoStart1, recoStart2, recoEnd1, recoEnd2;
1371 recoStart1.SetXYZ(track1->
Start().X(),
1372 track1->
Start().Y(),
1373 track1->
Start().Z());
1375 recoEnd1.SetXYZ(track1->
Stop().X(),
1377 track1->
Stop().Z());
1381 TVector3 highest1, lowest1;
1382 if (recoStart1.Y() > recoEnd1.Y())
1384 highest1 = recoStart1;
1393 recoStart2.SetXYZ(track2->
Start().X(),
1394 track2->
Start().Y(),
1395 track2->
Start().Z());
1397 recoEnd2.SetXYZ(track2->
Stop().X(),
1399 track2->
Stop().Z());
1401 TVector3 highest2, lowest2;
1402 if (recoStart2.Y() > recoEnd2.Y())
1404 highest2 = recoStart2;
1413 TVector3 dirvec1, dirvec2, diffvec;
1415 dirvec1 = highest1-lowest1;
1416 dirvec2 = highest2 - lowest2;
1451 unsigned int sumcross = YZcrosscount + XZcrosscount;
1453 unsigned int sumnocross = noXZcross + noYZcross;
1455 if (sumcross < sumnocross || sumcross == 0)
1459 nocrosstracks.push_back(track1);
1462 if (sumcross > sumnocross)
1465 std::cout <<
"TRACK " <<
i <<
" INTERSECTS " << sumcross <<
" TIMES IN XZ OR YZ VIEW. " <<
std::endl;
1474 if (nocrosstracks.size() < 2)
continue;
1476 fnInt += nocrosstracks.size();
1478 for (
unsigned int i = 0;
i < nocrosstracks.size(); ++
i)
1488 z2north(unitvec, unitvec_newxz);
1492 nocrossveccol.push_back(
trkvec);
1496 std::cout <<
"............Now Calculating DeltaTheta of NoCrossTracks." <<
std::endl;
1500 for (
unsigned int i=0;
i<nocrosstracks.size(); ++
i)
1508 z2north(unitvec, unitvec_newxz);
1513 std::vector<double> sigmacol;
1516 for (
unsigned int j =
i+1;
j < (
trkcol.size() - 1);
j++) {
1518 deltatheta = (
acos(scalarproduct))*180/
M_PI;
1539 std::vector<art::Ptr<rb::Track > > samesigncol;
1544 for (
unsigned int i = 0;
i < nocrosstracks.size(); ++
i)
1551 std::vector<art::Ptr<rb::Track > > tracks = nocrosstracks;
1554 std::vector<art::Ptr<rb::Track> > trkcollection = nocrosstracks;
1556 unsigned int oppsigntrks = 0;
1557 unsigned int samesigntrks = 0;
1562 if (oppsigntrks > samesigntrks) {
1564 std::cout <<
"Track " <<
i <<
" on slice " <<
nslices <<
"has opposite slope to " << oppsigntrks <<
" other tracks. " <<
std::endl;
1571 samesigncol.push_back(ftrack);
1580 if (samesigncol.size() < 2)
continue;
1582 fnOpp += samesigncol.size();
1584 for (
unsigned int i = 0;
i < samesigncol.size(); ++
i)
1594 z2north(unitvec, unitvec_newxz);
1598 samesignveccol.push_back(
trkvec);
1602 std::cout <<
"............Now Calculating DeltaTheta of SameSignTracks." <<
std::endl;
1604 for (
unsigned int i=0;
i<samesigncol.size(); ++
i)
1612 z2north(unitvec, unitvec_newxz);
1617 for (
unsigned int j =
i+1;
j < (
trkcol.size() - 1);
j++) {
1619 deltatheta = (
acos(scalarproduct))*180/
M_PI;
1624 minlengthtracks.push_back(track);
1628 if (minlengthtracks.size() < 2)
continue;
1629 fnLen += minlengthtracks.size();
1631 for (
unsigned int i = 0;
i < minlengthtracks.size(); ++
i)
1641 z2north(unitvec, unitvec_newxz);
1645 minlengthveccol.push_back(
trkvec);
1649 std::cout <<
"............Now Calculating DeltaTheta of MinLengthTracks." <<
std::endl;
1652 for (
unsigned int i=0;
i < minlengthtracks.size();
i++)
1660 z2north(unitvec, unitvec_newxz);
1663 trkcol = minlengthveccol;
1665 for (
unsigned int j =
i+1;
j < (minlengthtracks.size() - 1);
j++) {
1667 deltatheta = (
acos(scalarproduct))*180/
M_PI;
1675 std::cout <<
".......... Now Calculating Closest Distance Between Tracks" <<
std::endl;
1677 std::vector< double > DOCAmincol;
1679 for (
unsigned int i = 0 ;
i < minlengthtracks.size() ; ++
i)
1683 std::vector<double> DOCAcol;
1684 for (
unsigned int j =
i + 1;
j < (minlengthtracks.size() - 1);
j++)
1688 TVector3 recoStart1, recoStart2, recoEnd1, recoEnd2;
1690 recoStart1.SetXYZ(track1->
Start().X(),
1691 track1->
Start().Y(),
1692 track1->
Start().Z());
1694 recoEnd1.SetXYZ(track1->
Stop().X(),
1696 track1->
Stop().Z());
1698 TVector3 highest1, lowest1;
1699 if (recoStart1.Y() > recoEnd1.Y())
1701 highest1 = recoStart1;
1710 recoStart2.SetXYZ(track2->
Start().X(),
1711 track2->
Start().Y(),
1712 track2->
Start().Z());
1714 recoEnd2.SetXYZ(track2->
Stop().X(),
1716 track2->
Stop().Z());
1718 TVector3 highest2, lowest2;
1719 if (recoStart2.Y() > recoEnd2.Y())
1721 highest2 = recoStart2;
1730 TVector3 dirvec1, dirvec2, diffvec;
1732 dirvec1 = highest1-lowest1;
1733 dirvec2 = highest2 - lowest2;
1735 diffvec = highest1 - highest2;
1739 double a = dirvec1.X()*dirvec1.X() + dirvec1.Y()*dirvec1.Y() + dirvec1.Z()*dirvec1.Z();
1740 double b = dirvec1.X()*dirvec2.X() + dirvec1.Y()*dirvec2.Y() + dirvec1.Z()*dirvec2.Z();
1742 double c = dirvec2.X()*dirvec2.X() + dirvec2.Y()*dirvec2.Y() + dirvec2.Z()*dirvec2.Z();
1744 double d = dirvec1.X()*diffvec.X() + dirvec1.Y()*diffvec.Y() + dirvec1.Z()*diffvec.Z();
1746 double e = dirvec2.X()*diffvec.X() + dirvec2.Y()*diffvec.Y() + dirvec2.Z()*diffvec.Z();
1748 double D = a*c - b*
b;
1750 double sc, sN, sD = D;
1751 double tc, tN, tD = D;
1797 else if ((-d + b) >
a)
1805 sc = (TMath::Abs(sN) <
SMALL_NUM ? 0.0 : sN / sD);
1806 tc = (TMath::Abs(tN) <
SMALL_NUM ? 0.0 : tN / tD);
1809 TVector3 dP = diffvec + (sc * dirvec1 ) - (tc * dirvec2);
1811 double closestdist = dP.Mag();
1814 DOCAcol.push_back(closestdist);
1824 double DOCAmin = 10000;
1825 for (
unsigned int k = 0; k < DOCAcol.size(); k++)
1827 if (DOCAcol[
i] < DOCAmin)
1828 DOCAmin = DOCAcol[
i];
1831 DOCAmincol.push_back(DOCAmin);
1833 double DOCAmeter = DOCAmin/100.0;
1834 fDOCA->Fill(DOCAmeter);
1848 std::vector<double> azimuthcol;
1849 std::vector<double> altcol;
1850 std::vector<double> sigmacol;
1852 for (
unsigned int i = 0;
i < minlengthtracks.size();
i++){
1860 z2north(unitvec, unitvec_newxz);
1865 trkcol = minlengthveccol;
1868 for (
unsigned int j =
i+1;
j < (
trkcol.size() - 1);
j++) {
1870 deltatheta = (
acos(scalarproduct))*180/
M_PI;
1871 double deltaangle = (deltatheta*deltatheta)/2;
1879 track->
Start().Z());
1885 TVector3 highest, lowest;
1897 TVector3
length = highest - lowest;
1909 azimuthcol.push_back(azimuth);
1910 altcol.push_back(alt);
1926 std::cout <<
"...............Finished Processing SLICE " <<
nslices <<
" with " << minlengthtracks.size() <<
" passed tracks in event " << e.
event() <<
std::endl;
1936 std::cout <<
"************************************************************************" <<
std::endl;
1967 mf::LogError(
"AirKalmanAna") <<
"No Cosmic Exposure Info";
1993 srt_[
"cosmicrate"]=cosmicrate;
2003 <<
" cosmic tracks selected in this job.\n" 2007 <<
" according to CosmicExposure.\n" 2009 <<
"Note, livetime is for FULL SUBRUN, " 2011 <<
"so rate is wrong if haven't run full subrun.\n" 2013 <<
"The cosmic rate is therefore " 2015 << cosmicrate <<
" Hz \n";
AirKalmanAna(fhicl::ParameterSet const &p)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
void Unitize(art::Ptr< rb::Track > &track, TVector3 &unitvec)
SubRunNumber_t subRun() const
back track the reconstruction to the simulation
art::ServiceHandle< geo::LiveGeometry > livegeom
std::map< std::string, TH1 * > h_
constexpr std::uint32_t timeLow() const
virtual void beginSubRun(const art::SubRun &sr)
std::vector< double > anglecol
std::vector< double > pramcol
std::vector< unsigned > BadEventcol
double DistToEdge(double *point, double detHalfWidth, double detHalfHeight, double detLength)
Find the distance from the given point to the edge of the detector.
unsigned n3Dtrackspassedinevent
void AvgTrkVec(std::vector< TVector3 > &trkcol, TVector3 &avgtrk)
unsigned HighEnergyThresh
constexpr std::uint32_t timeHigh() const
void AngleBetween(TVector3 &avgtrk, TVector3 &trkvec, double &sigma)
unsigned n3Dtrackspassedinslice
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
void CheckAssociations(art::Ptr< rb::Cluster > slice, std::vector< art::Ptr< rb::Track > > tracks, bool filtered)
std::map< std::string, double > et_
void analyze(art::Event const &e) override
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
void AvgValue(std::vector< double > &pramcol, double &avgpram)
bool DoTracksIntersectXZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const
virtual double TotalLength() const
Length (cm) of all the track segments.
double TotalADC() const
Sum of the ADC of all the contained hits.
void OppSignTrkCount(art::Ptr< rb::Track > &ftrack, std::vector< art::Ptr< rb::Track > > &tracks, unsigned int &oppsigntrks, unsigned int &samesigntrks, unsigned int &i)
std::map< std::string, double > srt_
TH1D * fTracksPerSliceFinal
unsigned n3Dtracksinslice
void push_back(Ptr< U > const &p)
unsigned n3Dtracksinevent
void StdDev(std::vector< double > &anglecol, double &stddev)
TH2D * fAzimuthMultiplicity
Collect Geo headers and supply basic geometry functions.
void AngleBetweentrk(std::vector< TVector3 > &trkcol, TVector3 &trkvec, double &sigma, std::vector< double > &sigmacol)
art::ServiceHandle< geo::Geometry > geom
EventNumber_t event() const
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
void GetThetaXYZ(art::Ptr< rb::Track > &track, double &thetaX, double &thetaY, double &thetaZ)
std::vector< TVector3 > trkcol
T * make(ARGS...args) const
double DetHalfWidth() const
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.
bool DoTracksIntersectYZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const
void z2north(std::vector< TVector3 > track_i, std::vector< TVector3 > track_o)
unsigned n3Dtracksdeltathetapassedinevent
#define LOG_VERBATIM(category)
TH1D * fDeltaAnglePerSlice
unsigned n3Dtracksdeltathetapassedinslice
bool IsNoise() const
Is the noise flag set?
TVector3 Stop() const
Position of the final trajectory point.
std::vector< unsigned > Geventcol
bool CheckContainer(std::vector< art::Ptr< rb::Track > > &trkcollection, art::Ptr< rb::Track > &ftrack) const
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const