20 #include "ifdh_art/IFBeamService/IFBeam_service.h" 21 #include "art_root_io/TFileService.h" 39 #include "TPolyMarker.h" 52 typedef std::numeric_limits< double >
dbl;
95 void MaskGlitches( std::vector<short> &, std::array<short,192> & );
111 std::unique_ptr<ifbeam_ns::BeamFolder>&);
299 std::unique_ptr<ifbeam_ns::BeamFolder>
bfp;
343 L1(p.
get<double>(
"L1")),
344 L2(p.
get<double>(
"L2")),
345 L3(p.
get<double>(
"L3")),
416 produces<std::vector<beam::ProtoDUNEBeamEvent>>();
422 std::vector< std::pair<std::string, std::string> > tempTypes = p.
get<std::vector< std::pair<std::string, std::string> >>(
"DeviceTypes");
423 fDeviceTypes = std::map<std::string, std::string>(tempTypes.begin(), tempTypes.end() );
426 std::vector< std::pair<std::string, double> > tempFiberDims = p.
get< std::vector<std::pair<std::string,double> > >(
"Dimension");
427 fFiberDimension = std::map<std::string, double>(tempFiberDims.begin(), tempFiberDims.end());
441 std::vector<double> theResult;
443 MF_LOG_INFO(
"BeamEvent") <<
"Trying to grab from folder: " << name <<
"\n";
444 MF_LOG_INFO(
"BeamEvent") <<
"At Time: " << time <<
"\n";
447 theResult = the_folder->GetNamedVector(time, name);
450 MF_LOG_INFO(
"BeamEvent") <<
"Successfully fetched " << time <<
"\n";
466 MF_LOG_INFO(
"BeamEvent") <<
"Getting Raw Decoder Info" <<
"\n";
474 uint64_t high = RDTS.GetTimeStamp_High();
475 uint64_t low = RDTS.GetTimeStamp_Low();
477 uint64_t joined = (high | low);
483 MF_LOG_INFO(
"BeamEvent") <<
"High: " << RDTS.GetTimeStamp_High() <<
"\n";
484 MF_LOG_INFO(
"BeamEvent") <<
"Low: " << RDTS.GetTimeStamp_Low() <<
"\n";
485 MF_LOG_INFO(
"BeamEvent") <<
"Raw Decoder Timestamp: " << joined <<
"\n";
491 long long RDTSTickSec = (
RDTSTime * 2) / (
int)(TMath::Power(10,8));
492 RDTSTickSec = RDTSTickSec * (
int)(TMath::Power(10,8)) / 2;
493 long long RDTSTickNano =
RDTSTime - RDTSTickSec;
498 RDTSSec = 20.e-9 * RDTSTickSec;
507 auto const PDTStampHandle = e.
getValidHandle< std::vector< dune::ProtoDUNETimeStamp > > (
"timingrawdecoder:daq");
508 std::vector< dune::ProtoDUNETimeStamp > PDTStampVec(*PDTStampHandle);
509 auto PDTStamp = PDTStampVec[0];
511 UInt_t ver = PDTStamp.getVersion();
512 double RunStart = 2.e-8*PDTStamp.getLastRunStart();
513 SpillStart = 2.e-8*PDTStamp.getLastSpillStart();
514 SpillEnd = 2.e-8*PDTStamp.getLastSpillEnd();
523 MF_LOG_INFO(
"BeamEvent") <<
"Version: " << ver <<
"\n";
524 MF_LOG_INFO(
"BeamEvent") <<
"Run: " << RunStart <<
"\n";
536 std::vector<double> acqStamp =
FetchAndReport(time,
"dip/acc/NORTH/NP04/POW/MBPL022699:acqStamp[]",
bfp);
538 if( acqStamp[0] < 300000000.0 ){
540 MF_LOG_INFO(
"BeamEvent") <<
"Warning: MBPL Spill Start is low " << acqStamp[0]
541 <<
"\nWill need to time in with S11\n";
557 MF_LOG_WARNING(
"BeamEvent") <<
"Could not get Spill time to time in\n";
565 MF_LOG_INFO(
"BeamEvent") <<
"Matching S11 To Gen" <<
"\n";
567 for(
size_t iT = 0; iT <
beamspill->GetNT0(); ++iT){
568 double GenTrigSec =
beamspill->GetT0(iT).first;
569 double GenTrigNano =
beamspill->GetT0(iT).second;
571 double diff = GenTrigSec -
s11Sec;
572 diff += 1.e-9*(GenTrigNano -
s11Nano);
577 MF_LOG_INFO(
"BeamEvent") <<
"Found matching S11 and GenTrig!" <<
"\n";
578 MF_LOG_INFO(
"BeamEvent") <<
"diff: " << diff <<
"\n";
589 MF_LOG_INFO(
"BeamEvent") <<
"Could not find matching time " <<
"\n";
598 MF_LOG_INFO(
"BeamEvent") <<
"Matching in time between Beamline and TPC!!!" <<
"\n";
599 for(
size_t iT = 0; iT <
beamspill->GetNT0(); ++iT){
601 double GenTrigSec =
beamspill->GetT0(iT).first;
602 double GenTrigNano =
beamspill->GetT0(iT).second;
606 long long RDTSTickSec = (
RDTSTime * 2) / (
int)(TMath::Power(10,8));
607 RDTSTickSec = RDTSTickSec * (
int)(TMath::Power(10,8)) / 2;
608 long long RDTSTickNano =
RDTSTime - RDTSTickSec;
614 double diffSec = RDTSTimeSec - GenTrigSec -
SpillOffset;
615 double diffNano = 1.e-09*(RDTSTimeNano - GenTrigNano);
618 MF_LOG_INFO(
"BeamEvent") <<
"RDTSTimeSec - GenTrigSec " << RDTSTimeSec - GenTrigSec <<
"\n";
619 MF_LOG_INFO(
"BeamEvent") <<
"diff: " << diffSec <<
" " << diffNano <<
"\n";
622 double diff = diffSec + diffNano;
631 MF_LOG_INFO(
"BeamEvent") <<
"FOUND MATCHING TIME!!!" <<
"\n";
632 MF_LOG_INFO(
"BeamEvent") <<
"diff: " << diff <<
"\n";
640 MF_LOG_INFO(
"BeamEvent") <<
"Could not find matching time " <<
"\n";
647 MF_LOG_INFO(
"BeamEvent") <<
"art Event is unmatched to Beam Spill " <<
"\n";
651 size_t activeTrigger =
beamspill->GetActiveTrigger();
656 MF_LOG_INFO(
"BeamEvent") <<
"Setting beam event for matched event" <<
"\n";
662 for(
size_t i = 0; i <
fDevices.size(); ++i){
687 MF_LOG_INFO(
"BeamEvent") <<
"Finished adding info to beamevt " <<
"\n";
767 uint64_t fetch_time = uint64_t(
RDTSTime * 2e-8 );
768 uint64_t fetch_time_down = uint64_t(
RDTSTime * 2e-8 );
790 MF_LOG_INFO(
"BeamEvent") <<
"New spill or forced new fetch. Getting new beamspill info" <<
"\n";
799 MF_LOG_INFO(
"BeamEvent") <<
"fetch_time: " << fetch_time <<
"\n";
808 MF_LOG_INFO(
"BeamEvent") <<
"xcet fetch_time: " << fetch_time <<
"\n";
869 MF_LOG_INFO(
"BeamEvent") <<
"First Event: Priming cache\n";
941 MF_LOG_INFO(
"BeamEvent") <<
"Same spill. Reusing beamspill info" <<
"\n";
962 TimeIn(e, fetch_time_down);
1057 std::unique_ptr<std::vector<beam::ProtoDUNEBeamEvent> > beamData(
new std::vector<beam::ProtoDUNEBeamEvent>);
1098 std::vector<std::string> monitors;
1101 for(
size_t id = 0;
id <
fDevices.size(); ++id){
1104 monitors.push_back(name);
1114 MF_LOG_INFO(
"BeamEvent") <<
"Getting S11 Info " <<
"\n";
1116 std::vector<double> coarseS11 =
FetchAndReport(time,
"dip/acc/NORTH/NP04/BI/XBTF/S11:coarse[]",
bfp);
1117 std::vector<double> fracS11 =
FetchAndReport(time,
"dip/acc/NORTH/NP04/BI/XBTF/S11:frac[]",
bfp);
1118 std::vector<double> secondsS11 =
FetchAndReport(time,
"dip/acc/NORTH/NP04/BI/XBTF/S11:seconds[]",
bfp);
1119 std::vector<double> timestampCountS11 =
FetchAndReport(time,
"dip/acc/NORTH/NP04/BI/XBTF/S11:timestampCount",
bfp);
1121 int s11Count = (
int)timestampCountS11[0];
1124 MF_LOG_INFO(
"BeamEvent") <<
"Found " << s11Count <<
" returned signals" <<
"\n";
1128 long long RDTSTickSec = (
RDTSTime * 2) / (
int)(TMath::Power(10,8));
1129 RDTSTickSec = RDTSTickSec * (
int)(TMath::Power(10,8)) / 2;
1130 long long RDTSTickNano =
RDTSTime - RDTSTickSec;
1137 for(
int i = 0; i < s11Count; ++i){
1138 double nano = 8.*coarseS11[i] + fracS11[i]/512.;
1140 double diffSec = secondsS11[2*i + 1] - RDTSTimeSec -
fOffsetTAI;
1144 MF_LOG_INFO(
"BeamEvent") << i <<
" diffSec " << diffSec <<
"\n";
1145 MF_LOG_INFO(
"BeamEvent") << i <<
" diffNano " << diffNano <<
"\n";
1148 double diff = diffSec + 1.e-9*diffNano;
1152 MF_LOG_INFO(
"BeamEvent") <<
"FOUND Match between S11 and RDTS" <<
"\n";
1167 std::vector<double> coarseGeneralTrigger;
1168 std::vector<double> fracGeneralTrigger;
1169 std::vector<double> acqStampGeneralTrigger;
1170 std::vector<double> secondsGeneralTrigger;
1171 std::vector<double> timestampCountGeneralTrigger;
1175 coarseGeneralTrigger =
FetchAndReport(time,
"dip/acc/NORTH/NP04/BI/XBTF/GeneralTrigger:coarse[]",
bfp);
1176 fracGeneralTrigger =
FetchAndReport(time,
"dip/acc/NORTH/NP04/BI/XBTF/GeneralTrigger:frac[]",
bfp);
1177 acqStampGeneralTrigger =
FetchAndReport(time,
"dip/acc/NORTH/NP04/BI/XBTF/GeneralTrigger:acqStamp[]",
bfp);
1178 secondsGeneralTrigger =
FetchAndReport(time,
"dip/acc/NORTH/NP04/BI/XBTF/GeneralTrigger:seconds[]",
bfp);
1179 timestampCountGeneralTrigger =
FetchAndReport(time,
"dip/acc/NORTH/NP04/BI/XBTF/GeneralTrigger:timestampCount",
bfp);
1182 MF_LOG_INFO(
"BeamEvent") <<
"timestampCounts: " << timestampCountGeneralTrigger[0] <<
"\n";
1186 uint64_t low = (uint64_t)acqStampGeneralTrigger[1];
1187 uint64_t high = (uint64_t)acqStampGeneralTrigger[0];
1189 acqTime = ( high | low ) / 1000000000.;
1193 MF_LOG_WARNING(
"BeamEvent") <<
"Could not get GeneralTrigger information!!" <<
"\n";
1198 std::vector<double> coarseTOF1A;
1199 std::vector<double> fracTOF1A;
1200 std::vector<double> secondsTOF1A;
1202 std::vector<double> coarseTOF1B;
1203 std::vector<double> fracTOF1B;
1204 std::vector<double> secondsTOF1B;
1206 std::vector<double> coarseTOF2A;
1207 std::vector<double> fracTOF2A;
1208 std::vector<double> secondsTOF2A;
1210 std::vector<double> coarseTOF2B;
1211 std::vector<double> fracTOF2B;
1212 std::vector<double> secondsTOF2B;
1234 MF_LOG_WARNING(
"BeamEvent") <<
"Could not get TOF information!!" <<
"\n";
1238 std::vector<std::pair<double,double>> unorderedGenTrigTime;
1239 std::vector<std::pair<double,double>> unorderedTOF1ATime;
1240 std::vector<std::pair<double,double>> unorderedTOF1BTime;
1241 std::vector<std::pair<double,double>> unorderedTOF2ATime;
1242 std::vector<std::pair<double,double>> unorderedTOF2BTime;
1245 for(
size_t i = 0; i < timestampCountGeneralTrigger[0]; ++i){
1248 MF_LOG_INFO(
"BeamEvent") << i <<
" " << secondsGeneralTrigger[2*i + 1] <<
" " << 8.*coarseGeneralTrigger[i] + fracGeneralTrigger[i]/512. <<
"\n";
1249 MF_LOG_INFO(
"BeamEvent") <<
"\t" <<
std::setw(15) << secondsGeneralTrigger[2*i + 1] + 1.e-9*(8.*coarseGeneralTrigger[i] + fracGeneralTrigger[i]/512.) <<
"\n";
1269 for(
size_t i = 0; i < timestampCountGeneralTrigger[0]; ++i){
1272 MF_LOG_INFO(
"BeamEvent") <<
"TOF1A " << i <<
" " << secondsTOF1A[2*i+1] <<
" " << 8.*coarseTOF1A[i] + fracTOF1A[i]/512. <<
"\n";
1286 for(
size_t i = 0; i < timestampCountGeneralTrigger[0]; ++i){
1289 MF_LOG_INFO(
"BeamEvent") <<
"TOF1B " << i <<
" " << secondsTOF1B[2*i+1] <<
" " << 8.*coarseTOF1B[i] + fracTOF1B[i]/512. <<
"\n";
1302 for(
size_t i = 0; i < timestampCountGeneralTrigger[0]; ++i){
1305 MF_LOG_INFO(
"BeamEvent") <<
"TOF2A " << i <<
" " << secondsTOF2A[2*i+1] <<
" " << 8.*coarseTOF2A[i] + fracTOF2A[i]/512. <<
"\n";
1317 for(
size_t j = 0; j < timestampCountGeneralTrigger[0]; ++j){
1318 diff2A.push_back( 1.e9*(unorderedTOF2ATime.back().first - unorderedGenTrigTime[j].first) + (unorderedTOF2ATime.back().second - unorderedGenTrigTime[j].second) );
1325 for(
size_t i = 0; i < timestampCountGeneralTrigger[0]; ++i){
1328 MF_LOG_INFO(
"BeamEvent") <<
"TOF2B " << i <<
" " << secondsTOF2B[2*i+1] <<
" " << 8.*coarseTOF2B[i] + fracTOF2B[i]/512. <<
"\n";
1341 for(
size_t j = 0; j < timestampCountGeneralTrigger[0]; ++j){
1342 diff2B.push_back( 1.e9*(unorderedTOF2BTime.back().first - unorderedGenTrigTime[j].first) + (unorderedTOF2BTime.back().second - unorderedGenTrigTime[j].second) );
1350 MF_LOG_INFO(
"BeamEvent") <<
"NGenTrigs: " << timestampCountGeneralTrigger[0] <<
" NTOF2s: " << unorderedTOF2ATime.size() + unorderedTOF2BTime.size() <<
"\n";
1352 for(
size_t iT = 0; iT < unorderedGenTrigTime.size(); ++iT){
1354 bool found_TOF =
false;
1356 double the_gen_sec = unorderedGenTrigTime[iT].first;
1357 double the_gen_ns = unorderedGenTrigTime[iT].second;
1363 std::vector< double > possibleTOF;
1364 std::vector< int > possibleTOFChan;
1365 std::vector< size_t > UpstreamTriggers;
1366 std::vector< size_t > DownstreamTriggers;
1371 MF_LOG_INFO(
"BeamEvent") <<
"Gen: " << the_gen_sec <<
" " << the_gen_ns <<
"\n";
1374 for(
size_t ip2A = 0; ip2A < unorderedTOF2ATime.size(); ++ip2A){
1375 double TOF2A_sec = unorderedTOF2ATime[ip2A].first;
1376 double TOF2A_ns = unorderedTOF2ATime[ip2A].second;
1377 double delta_2A = 1.e9*(the_gen_sec - TOF2A_sec) + the_gen_ns - TOF2A_ns;
1379 if( delta_2A < 0. )
break;
1384 MF_LOG_INFO(
"BeamEvent") <<
"Found match 2A to Gen" <<
"\n";
1387 for(
size_t ip1A = 0; ip1A < unorderedTOF1ATime.size(); ++ip1A){
1388 double TOF1A_sec = unorderedTOF1ATime[ip1A].first;
1389 double TOF1A_ns = unorderedTOF1ATime[ip1A].second;
1390 double delta = 1.e9*( TOF2A_sec - TOF1A_sec ) + TOF2A_ns - TOF1A_ns;
1392 if( delta < 0. )
break;
1395 MF_LOG_INFO(
"BeamEvent") <<
"Found match 1A to 2A " << delta <<
"\n";
1400 possibleTOF.push_back( delta );
1401 possibleTOFChan.push_back( channel );
1403 UpstreamTriggers.push_back( ip1A );
1404 DownstreamTriggers.push_back( ip2A );
1408 for(
size_t ip1B = 0; ip1B < unorderedTOF1BTime.size(); ++ip1B){
1409 double TOF1B_sec = unorderedTOF1BTime[ip1B].first;
1410 double TOF1B_ns = unorderedTOF1BTime[ip1B].second;
1411 double delta = 1.e9*( TOF2A_sec - TOF1B_sec ) + TOF2A_ns - TOF1B_ns;
1413 if( delta < 0. )
break;
1417 MF_LOG_INFO(
"BeamEvent") <<
"Found match 1B to 2A " << delta <<
"\n";
1422 possibleTOF.push_back( delta );
1423 possibleTOFChan.push_back( channel );
1425 UpstreamTriggers.push_back( ip1B );
1426 DownstreamTriggers.push_back( ip2A );
1433 for(
size_t ip2B = 0; ip2B < unorderedTOF2BTime.size(); ++ip2B){
1434 double TOF2B_sec = unorderedTOF2BTime[ip2B].first;
1435 double TOF2B_ns = unorderedTOF2BTime[ip2B].second;
1436 double delta_2B = 1.e9*(the_gen_sec - TOF2B_sec) + the_gen_ns - TOF2B_ns;
1439 if( delta_2B < 0. )
break;
1444 MF_LOG_INFO(
"BeamEvent") <<
"Found match 2B to Gen" <<
"\n";
1447 for(
size_t ip1A = 0; ip1A < unorderedTOF1ATime.size(); ++ip1A){
1448 double TOF1A_sec = unorderedTOF1ATime[ip1A].first;
1449 double TOF1A_ns = unorderedTOF1ATime[ip1A].second;
1450 double delta = 1.e9*( TOF2B_sec - TOF1A_sec ) + TOF2B_ns - TOF1A_ns;
1453 if( delta < 0. )
break;
1457 MF_LOG_INFO(
"BeamEvent") <<
"Found match 1A to 2B " << delta <<
"\n";
1462 possibleTOF.push_back( delta );
1463 possibleTOFChan.push_back( channel );
1465 UpstreamTriggers.push_back( ip1A );
1466 DownstreamTriggers.push_back( ip2B );
1470 for(
size_t ip1B = 0; ip1B < unorderedTOF1BTime.size(); ++ip1B){
1471 double TOF1B_sec = unorderedTOF1BTime[ip1B].first;
1472 double TOF1B_ns = unorderedTOF1BTime[ip1B].second;
1473 double delta = 1.e9*( TOF2B_sec - TOF1B_sec ) + TOF2B_ns - TOF1B_ns;
1476 if( delta < 0. )
break;
1480 MF_LOG_INFO(
"BeamEvent") <<
"Found match 1B to 2B " << delta <<
"\n";
1485 possibleTOF.push_back( delta );
1486 possibleTOFChan.push_back( channel );
1488 UpstreamTriggers.push_back( ip1B );
1489 DownstreamTriggers.push_back( ip2B );
1496 MF_LOG_INFO(
"BeamEvent") <<
"Found " << possibleTOF.size() <<
" matched TOFs" <<
"\n";
1502 MF_LOG_INFO(
"BeamEvent") <<
"No matching TOFs found. Placing dummy\n";
1504 possibleTOF.push_back( 0. );
1505 possibleTOFChan.push_back( -1 );
1506 UpstreamTriggers.push_back(0);
1507 DownstreamTriggers.push_back(0);
1511 beamspill->AddMultipleTOFs( possibleTOF );
1512 beamspill->AddMultipleTOFChans( possibleTOFChan );
1513 beamspill->AddUpstreamTriggers( UpstreamTriggers );
1514 beamspill->AddDownstreamTriggers( DownstreamTriggers );
1534 MF_LOG_WARNING(
"BeamEvent") <<
"Could not get Cerenkov 1 Pressure\n";
1549 MF_LOG_WARNING(
"BeamEvent") <<
"Could not get Cerenkov 2 Pressure\n";
1555 std::vector< double > XCET1_timestamps, XCET2_timestamps;
1556 std::vector< double > XCET1_seconds, XCET2_seconds;
1557 std::vector< double > XCET1_frac, XCET2_frac;
1558 std::vector< double > XCET1_coarse, XCET2_coarse;
1560 bool fetched_XCET1=
false;
1561 bool fetched_XCET2=
false;
1567 fetched_XCET1 =
true;
1570 for(
size_t i = 0; i < XCET1_seconds.size(); ++i ){
1571 std::cout << i <<
" " << XCET1_seconds[i] -
fOffsetTAI <<
" " << (8.*XCET1_coarse[i] + XCET1_frac[i] / 512.) <<
std::endl;
1578 fetched_XCET1 =
false;
1586 fetched_XCET2 =
true;
1589 for(
size_t i = 0; i < XCET2_seconds.size(); ++i ){
1590 std::cout << i <<
" " << XCET2_seconds[i] -
fOffsetTAI <<
" " << (8.*XCET2_coarse[i] + XCET2_frac[i] / 512.) <<
std::endl;
1597 fetched_XCET2 =
false;
1604 for(
size_t i = 0; i <
beamspill->GetNT0(); ++i ){
1612 if( !fetched_XCET1 ) status_1.
trigger = -1;
1617 for(
size_t ic1 = 0; ic1 < XCET1_seconds.size(); ++ic1 ){
1619 delta += (
beamspill->GetT0Nano(i) - (8.*XCET1_coarse[ic1] + XCET1_frac[ic1] / 512.) );
1623 if( fabs(delta) < 500. ){
1624 if(
fXCETDebug ) std::cout <<
"Found matching XCET1 trigger " << XCET1_seconds[ic1] -
fOffsetTAI <<
" " << (8.*XCET1_coarse[ic1] + XCET1_frac[ic1] / 512.) <<
" " << delta <<
std::endl;
1637 if( !fetched_XCET2 ) status_2.
trigger = -1;
1642 for(
size_t ic2 = 0; ic2 < XCET2_seconds.size(); ++ic2 ){
1645 delta += (
beamspill->GetT0Nano(i) - (8.*XCET2_coarse[ic2] + XCET2_frac[ic2] / 512.) );
1649 if( fabs(delta) < 500. ){
1650 if(
fXCETDebug ) std::cout <<
"Found matching XCET2 trigger " << XCET2_seconds[ic2] -
fOffsetTAI <<
" " << (8.*XCET2_coarse[ic2] + XCET2_frac[ic2] / 512.) <<
" " << delta <<
std::endl;
1666 int nxcet1 = 0, nxcet2 = 0;
1668 for(
size_t i = 0; i <
beamspill->GetNT0(); ++i ){
1669 if(
beamspill->GetCKov0Status(i) == 1 ) nxcet1++ ;
1670 if(
beamspill->GetCKov1Status(i) == 1 ) nxcet2++ ;
1675 if( fetched_XCET1 )
MF_LOG_INFO(
"BeamEvent") <<
"XCET1: " << XCET1_seconds.size() <<
" " << nxcet1 <<
std::endl;
1676 if( fetched_XCET2 )
MF_LOG_INFO(
"BeamEvent") <<
"XCET2: " << XCET2_seconds.size() <<
" " << nxcet2 <<
std::endl;
1687 std::vector<double>
counts;
1698 MF_LOG_INFO(
"BeamEvent") <<
"Counts: " << counts[1] <<
"\n";
1700 std::vector<double>
data;
1733 std::vector<size_t> leftOvers;
1734 for(
size_t lo = 0; lo <
beamspill->GetNT0(); ++lo){
1735 leftOvers.push_back(lo);
1739 for(
size_t i = 0; i < counts[1]; ++i){
1740 for(
int j = 0; j < 10; ++j){
1742 double theData = data[20*i + (2*j + 1)];
1765 if( delta < -1.e-7 && delta > -10.
e-7 ){
1767 beamspill->ReplaceFBMTrigger(name, fbm, *ip);
1768 leftOvers.erase(ip);
1775 if( leftOvers.size() ){
1776 MF_LOG_WARNING(
"BeamEvent") <<
"Warning! Could not match to Good Particles: " <<
"\n";
1777 for(
size_t ip = 0; ip < leftOvers.size(); ++ip){
1784 for(
size_t i = 0; i <
beamspill->GetNFBMTriggers(name); ++i){
1791 std::cout <<
"Fixed " << name <<
std::endl;
1792 for(
size_t i = 0; i <
beamspill->GetNFBMTriggers(name); ++i ){
1793 for(
size_t j = 0; j < 192; ++j ){
1794 if(
beamspill->GetFBM( name, i ).glitch_mask[j] ){
1795 std::cout << i <<
" " << j <<
std::endl;
1821 fFullMomentum = tfs->make<TH1F>(
"FullMomentum",
"", 150, 0, 15.);
1822 fCutMomentum = tfs->make<TH1F>(
"CutMomentum",
"", 150, 0, 15.);
1826 fOutTree = tfs->make<TTree>(
"tree",
"lines");
1894 MF_LOG_INFO(
"BeamEvent") <<
"%%%%%%%%%% Got beam folder %%%%%%%%%%\n";
1906 MF_LOG_INFO(
"BeamEvent") <<
"%%%%%%%%%% Got beam folder %%%%%%%%%%\n";
1922 uint64_t low64 = (uint64_t)low;
1924 uint64_t high64 = (uint64_t)high;
1926 high64 = high64 << 32;
1927 uint64_t joined = high64 | low64;
1935 TVector3 old(x,y,z);
1945 TVector3
result(newX/10., newY/10., newZ/10.);
1965 MF_LOG_INFO(
"BeamEvent") <<
"Making Track for time: " <<
beamspill->GetT0Sec(theTrigger) <<
" " <<
beamspill->GetT0Nano(theTrigger) <<
"\n";
1973 for(
size_t i = 0; i < firstUpstreamFibers.size(); ++i){
1974 MF_LOG_INFO(
"BeamEvent") << firstUpstreamFibers[i] <<
" ";
1979 for(
size_t i = 0; i < secondUpstreamFibers.size(); ++i){
1980 MF_LOG_INFO(
"BeamEvent") << secondUpstreamFibers[i] <<
" ";
1987 MaskGlitches( firstUpstreamFibers, firstUpstreamGlitches );
1988 MaskGlitches( secondUpstreamFibers, secondUpstreamGlitches );
1999 for(
size_t i = 0; i < firstDownstreamFibers.size(); ++i){
2000 MF_LOG_INFO(
"BeamEvent") << firstDownstreamFibers[i] <<
" ";
2005 for(
size_t i = 0; i < secondDownstreamFibers.size(); ++i){
2006 MF_LOG_INFO(
"BeamEvent") << secondDownstreamFibers[i] <<
" ";
2012 MaskGlitches( firstDownstreamFibers, firstDownstreamGlitches );
2013 MaskGlitches( secondDownstreamFibers, secondDownstreamGlitches );
2018 if( (firstUpstreamFibers.size() < 1) || (secondUpstreamFibers.size() < 1) || (firstDownstreamFibers.size() < 1) || (secondDownstreamFibers.size() < 1) ){
2020 MF_LOG_INFO(
"BeamEvent") <<
"Warning, at least one empty Beam Profiler. Not making track" <<
"\n";
2026 std::vector< std::pair<size_t, size_t> > upstreamPairedFibers;
2027 std::vector< std::pair<size_t, size_t> > downstreamPairedFibers;
2033 std::vector< TVector3 > upstreamPositions;
2034 std::vector< TVector3 > downstreamPositions;
2040 for(
size_t iF1 = 0; iF1 < firstUpstreamFibers.size(); ++iF1){
2042 size_t firstFiber = firstUpstreamFibers[iF1];
2044 for(
size_t iF2 = 0; iF2 < secondUpstreamFibers.size(); ++iF2){
2045 size_t secondFiber = secondUpstreamFibers[iF2];
2048 MF_LOG_INFO(
"BeamEvent") <<
"Paired: " << firstFiber <<
" " << secondFiber <<
"\n";
2050 upstreamPairedFibers.push_back(std::make_pair(firstFiber, secondFiber));
2052 if (iF2 < secondUpstreamFibers.size() - 1){
2053 if (secondUpstreamFibers[iF2] == (secondUpstreamFibers[iF2 + 1] - 1)) ++iF2;
2057 if (iF1 < firstUpstreamFibers.size() - 1){
2058 if (firstUpstreamFibers[iF1] == (firstUpstreamFibers[iF1 + 1] - 1)) ++iF1;
2062 for(
size_t iF = 0; iF < upstreamPairedFibers.size(); ++iF){
2064 std::pair<size_t,size_t> thePair = upstreamPairedFibers.at(iF);
2066 if(firstUpstreamType ==
"horiz" && secondUpstreamType ==
"vert"){
2071 upstreamPositions.push_back( posInDet );
2074 MF_LOG_INFO(
"BeamEvent") <<
"normal " << xPos <<
" " << yPos <<
"\n";
2075 MF_LOG_INFO(
"BeamEvent") << posInDet.X() <<
" " << posInDet.Y() <<
" " << posInDet.Z() <<
"\n";
2079 else if(firstUpstreamType ==
"vert" && secondUpstreamType ==
"horiz"){
2084 upstreamPositions.push_back( posInDet );
2087 MF_LOG_INFO(
"BeamEvent") <<
"normal " << xPos <<
" " << yPos <<
"\n";
2088 MF_LOG_INFO(
"BeamEvent") << posInDet.X() <<
" " << posInDet.Y() <<
" " << posInDet.Z() <<
"\n";
2098 for(
size_t iF1 = 0; iF1 < firstDownstreamFibers.size(); ++iF1){
2100 size_t firstFiber = firstDownstreamFibers[iF1];
2102 for(
size_t iF2 = 0; iF2 < secondDownstreamFibers.size(); ++iF2){
2103 size_t secondFiber = secondDownstreamFibers[iF2];
2106 MF_LOG_INFO(
"BeamEvent") <<
"Paired: " << firstFiber <<
" " << secondFiber <<
"\n";
2107 downstreamPairedFibers.push_back(std::make_pair(firstFiber, secondFiber));
2109 if (iF2 < secondDownstreamFibers.size() - 1){
2110 if (secondDownstreamFibers[iF2] == (secondDownstreamFibers[iF2 + 1] - 1)) ++iF2;
2114 if (iF1 < firstDownstreamFibers.size() - 1){
2115 if (firstDownstreamFibers[iF1] == (firstDownstreamFibers[iF1 + 1] - 1)) ++iF1;
2119 for(
size_t iF = 0; iF < downstreamPairedFibers.size(); ++iF){
2121 std::pair<size_t,size_t> thePair = downstreamPairedFibers.at(iF);
2123 if(firstDownstreamType ==
"horiz" && secondDownstreamType ==
"vert"){
2128 downstreamPositions.push_back( posInDet );
2131 MF_LOG_INFO(
"BeamEvent") <<
"normal " << xPos <<
" " << yPos <<
"\n";
2132 MF_LOG_INFO(
"BeamEvent") << posInDet.X() <<
" " << posInDet.Y() <<
" " << posInDet.Z() <<
"\n";
2135 else if(firstDownstreamType ==
"vert" && secondDownstreamType ==
"horiz"){
2140 downstreamPositions.push_back( posInDet );
2143 MF_LOG_INFO(
"BeamEvent") <<
"normal " << xPos <<
" " << yPos <<
"\n";
2144 MF_LOG_INFO(
"BeamEvent") << posInDet.X() <<
" " << posInDet.Y() <<
" " << posInDet.Z() <<
"\n";
2152 for(
size_t iU = 0; iU < upstreamPositions.size(); ++iU){
2153 for(
size_t iD = 0; iD < downstreamPositions.size(); ++iD){
2154 std::vector<TVector3> thePoints;
2155 thePoints.push_back(upstreamPositions.at(iU));
2156 thePoints.push_back(downstreamPositions.at(iD));
2159 thePoints.push_back(
ProjectToTPC(thePoints[0],thePoints[1]) );
2160 std::vector<TVector3> theMomenta;
2162 theMomenta.push_back( ( downstreamPositions.at(iD) - upstreamPositions.at(iU) ).Unit() );
2163 theMomenta.push_back( ( downstreamPositions.at(iD) - upstreamPositions.at(iU) ).Unit() );
2164 theMomenta.push_back( ( downstreamPositions.at(iD) - upstreamPositions.at(iU) ).Unit() );
2180 MF_LOG_INFO(
"BeamEvent") <<
"Doing momentum spectrometry for trigger " <<
beamspill->GetT0Sec(theTrigger) <<
" " <<
beamspill->GetT0Nano(theTrigger) <<
"\n";
2189 std::vector<short> BPROF1Fibers;
2192 if (firstBPROF1Type ==
"horiz" && secondBPROF1Type ==
"vert"){
2196 else if(secondBPROF1Type ==
"horiz" && firstBPROF1Type ==
"vert"){
2201 MF_LOG_WARNING(
"BeamEvent") <<
"Error: type is not correct" <<
"\n";
2206 MF_LOG_INFO(
"BeamEvent") << BPROF1Name <<
" has " << BPROF1Fibers.size() <<
" active fibers" <<
"\n";
2207 std::cout.precision(20);
2208 std::cout <<
"Data: \n" 2209 <<
"\t" <<
beamspill->GetFBM( BPROF1Name, theTrigger ).fiberData[0] <<
"\n" 2210 <<
"\t" <<
beamspill->GetFBM( BPROF1Name, theTrigger ).fiberData[1] <<
"\n" 2211 <<
"\t" <<
beamspill->GetFBM( BPROF1Name, theTrigger ).fiberData[2] <<
"\n" 2212 <<
"\t" <<
beamspill->GetFBM( BPROF1Name, theTrigger ).fiberData[3] <<
"\n" 2213 <<
"\t" <<
beamspill->GetFBM( BPROF1Name, theTrigger ).fiberData[4] <<
"\n" 2216 for(
size_t i = 0; i < BPROF1Fibers.size(); ++i){
2217 MF_LOG_INFO(
"BeamEvent") << BPROF1Fibers[i] <<
" ";
2224 std::array<short,192> BPROF1Glitches =
beamspill->GetFBM( BPROF1Name, theTrigger ).glitch_mask;
2225 if(
fPrintDebug ) std::cout <<
"Got " << BPROF1Fibers.size() <<
" Fibers before Fix" <<
std::endl;
2231 if( (BPROF1Fibers.size() < 1) ){
2233 MF_LOG_INFO(
"BeamEvent") <<
"Warning, at least one empty Beam Profiler. Not checking momentum" <<
"\n";
2241 std::vector<short> BPROF2Fibers =
beamspill->GetActiveFibers(
BPROF2, theTrigger);
2244 std::array<short,192> BPROF2Glitches =
beamspill->GetFBM(
BPROF2, theTrigger ).glitch_mask;
2245 if(
fPrintDebug ) std::cout <<
"Got " << BPROF2Fibers.size() <<
" Fibers before Fix" <<
std::endl;
2249 if( (BPROF2Fibers.size() < 1) ){
2251 MF_LOG_INFO(
"BeamEvent") <<
"Warning, at least one empty Beam Profiler. Not checking momentum" <<
"\n";
2257 for(
size_t i = 0; i < BPROF2Fibers.size(); ++i){
2258 MF_LOG_INFO(
"BeamEvent") << BPROF2Fibers[i] <<
" ";
2266 std::vector<short> BPROF3Fibers =
beamspill->GetActiveFibers(
BPROF3, theTrigger);
2269 std::array<short,192> BPROF3Glitches =
beamspill->GetFBM(
BPROF3, theTrigger ).glitch_mask;
2270 if(
fPrintDebug ) std::cout <<
"Got " << BPROF3Fibers.size() <<
" Fibers before Fix" <<
std::endl;
2275 if( (BPROF3Fibers.size() < 1) ){
2277 MF_LOG_INFO(
"BeamEvent") <<
"Warning, at least one empty Beam Profiler. Not checking momentum" <<
"\n";
2283 for(
size_t i = 0; i < BPROF3Fibers.size(); ++i){
2284 MF_LOG_INFO(
"BeamEvent") << BPROF3Fibers[i] <<
" ";
2290 MF_LOG_INFO(
"BeamEvent") <<
"Getting all trio-wise hits" <<
"\n";
2291 MF_LOG_INFO(
"BeamEvent") <<
"N1,N2,N3 " << BPROF1Fibers.size()
2292 <<
" " << BPROF2Fibers.size()
2293 <<
" " << BPROF3Fibers.size() <<
"\n";
2296 for(
size_t i1 = 0; i1 < BPROF1Fibers.size(); ++i1){
2299 x1 = -1.*
GetPosition(BPROF1Name, BPROF1Fibers[i1])/1.E3;
2300 if (i1 < BPROF1Fibers.size() - 1){
2301 if (BPROF1Fibers[i1] == (BPROF1Fibers[i1 + 1] - 1)){
2307 for(
size_t i2 = 0; i2 < BPROF2Fibers.size(); ++i2){
2309 if (i2 < BPROF2Fibers.size() - 1){
2310 if (BPROF2Fibers[i2] == (BPROF2Fibers[i2 + 1] - 1)){
2316 for(
size_t i3 = 0; i3 < BPROF3Fibers.size(); ++i3){
2319 MF_LOG_INFO(
"BeamEvent") <<
"\t" << i1 <<
" " << i2 <<
" " << i3 <<
"\n";
2322 if (i3 < BPROF3Fibers.size() - 1){
2323 if (BPROF3Fibers[i3] == (BPROF3Fibers[i3 + 1] - 1)){
2338 double momentum_full = 299792458*LB/(1.E9 * acos(cosTheta_full));
2343 if (i3 < BPROF3Fibers.size() - 1){
2344 if (BPROF3Fibers[i3] == (BPROF3Fibers[i3 + 1] - 1)){
2352 if (i2 < BPROF2Fibers.size() - 1){
2353 if (BPROF2Fibers[i2] == (BPROF2Fibers[i2 + 1] - 1)){
2360 if (i1 < BPROF1Fibers.size() - 1){
2361 if (BPROF1Fibers[i1] == (BPROF1Fibers[i1 + 1] - 1)){
2375 double denomTerm1, denomTerm2, denom;
2376 denomTerm1 = sqrt(
L1*
L1 + (a - X1)*(a - X1) );
2378 + TMath::Power( ( (
L3 -
L2) ),2) );
2379 denom = denomTerm1 * denomTerm2;
2381 double cosTheta = numTerm/denom;
2386 TVector3 dR = (secondPoint - firstPoint);
2388 double deltaZ = -1.*secondPoint.Z();
2389 double deltaX = deltaZ * (dR.X() / dR.Z());
2390 double deltaY = deltaZ * (dR.Y() / dR.Z());
2392 TVector3 lastPoint = secondPoint + TVector3(deltaX, deltaY, deltaZ);
2397 if(fiberIdx > 192){
MF_LOG_WARNING(
"BeamEvent") <<
"Please select fiber in range [0,191]" <<
"\n";
return -1.;}
2401 double pos = size*(96 - fiberIdx) - size/2.;
2406 for(
short i = 0; i < 192; ++i ){
2408 fibers.erase( std::find( fibers.begin(), fibers.end(), i ) );
void GetSpillInfo(art::Event &)
std::string fXCETBundleName
double GetPosition(std::string, int)
void InitXBPFInfo(beam::ProtoDUNEBeamSpill *)
constexpr std::uint32_t timeLow() const
EventNumber_t event() const
const double & GetT0Nano() const
size_t GetNRecoBeamMomenta()
std::vector< std::string > fDevices
const int & GetTimingTrigger() const
void SetSpillStart(double theSpillStart)
void ClearRecoBeamMomenta()
std::vector< double > FetchAndReport(long long, std::string, std::unique_ptr< ifbeam_ns::BeamFolder > &)
Handle< PROD > getHandle(SelectorBase const &) const
void getS11Info(uint64_t)
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 > > SMatrixSym55
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
const double & GetCKov0Pressure() const
std::map< std::string, double > fFiberDimension
std::vector< double > diff2B
void SetTimingTrigger(int theTrigger)
EDProducer(fhicl::ParameterSet const &pset)
constexpr std::uint32_t timeHigh() const
std::string firstDownstreamName
std::unique_ptr< ifbeam_ns::BeamFolder > bfp_xcet
BeamEvent(fhicl::ParameterSet const &p)
const short & GetCKov1Status() const
TrackTrajectory::Flags_t Flags_t
double fDownstreamToGenTrig
void MaskGlitches(std::vector< short > &, std::array< short, 192 > &)
#define MF_LOG_ERROR(category)
double GetPairedPosition(std::string, size_t)
art::Handle< std::vector< raw::RDTimeStamp > > RDTimeStampHandle
double fFirstTrackingProfZ
double fUpstreamToDownstream
beam::ProtoDUNEBeamEvent prev_beamevt
void InitFBMs(std::vector< std::string >)
TVector3 ProjectToTPC(TVector3, TVector3)
beam::ProtoDUNEBeamEvent * beamevt
void parseXCETDB(uint64_t)
void BeamMonitorBasisVectors()
BeamEvent & operator=(BeamEvent const &)=delete
void SetBITrigger(int theTrigger)
void SetSpillOffset(double theSpillOffset)
std::unique_ptr< ifbeam_ns::BeamFolder > bfp
void SetFBMTrigger(std::string, FBM)
void MomentumSpec(size_t)
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
void SetCKov0(CKov theCKov)
double MomentumCosTheta(double, double, double)
#define DEFINE_ART_MODULE(klass)
void SetMagnetCurrent(double theMagnetCurrent)
std::vector< double > genTrigCoarses
A trajectory in space reconstructed from hits.
const double & GetTOF() const
const double & GetT0Sec() const
counts_as<> counts
Number of ADC counts, represented by signed short int.
art::ServiceHandle< ifbeam_ns::IFBeam > ifb
T get(std::string const &key) const
std::vector< double > diff2A
void SetRDTimestamp(long long theRDTimestamp)
double fCalibrationTolerance
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
uint64_t joinHighLow(double, double)
SubRunNumber_t subRun() const
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
void SetCTBTimestamp(long long theCTBTimestamp)
General LArSoft Utilities.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
void AddBeamTrack(recob::Track theTrack)
#define MF_LOG_INFO(category)
std::string firstUpstreamName
std::vector< double > current
beam::ProtoDUNEBeamSpill * beamspill
Q_EXPORT QTSManip setw(int w)
void AddRecoBeamMomentum(double theMomentum)
double fSecondTrackingProfZ
double fTimingCalibration
const size_t & GetActiveTrigger() const
std::unique_ptr< BeamFolder > getBeamFolder(std::string bundle_name, std::string url="", double time_width=1200.0)
std::map< std::string, std::string > fDeviceTypes
void SetDownstreamTriggers(std::vector< size_t > theContent)
std::string secondDownstreamName
Provides recob::Track data product.
void SetUpstreamTriggers(std::vector< size_t > theContent)
void TimeIn(art::Event &, uint64_t)
const short & GetCKov0Status() const
void SetT0(std::pair< double, double > theT0)
void parseGeneralXBPF(std::string, uint64_t, size_t)
void SetCKov1(CKov theCKov)
auto const & get(AssnsNode< L, R, D > const &r)
#define MF_LOG_WARNING(category)
beam::ProtoDUNEBeamSpill prev_beamspill
TVector3 ConvertProfCoordinates(double x, double y, double z, double zOffset)
const double & GetRecoBeamMomentum(size_t i) const
const double & GetMagnetCurrent() const
const bool & CheckIsMatched() const
void RotateMonitorVector(TVector3 &vec)
void SetTOFs(std::vector< double > theContent)
std::vector< double > genTrigFracs
const int & GetTOFChan() const
std::numeric_limits< double > dbl
void SetCalibrations(double TOFCalAA, double TOFCalBA, double TOFCalAB, double TOFCalBB)
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
std::vector< double > genTrigSecs
cet::coded_exception< error, detail::translate > exception
std::string secondUpstreamName
void SetTOFChans(std::vector< int > theContent)
void produce(art::Event &e) override
QTextStream & endl(QTextStream &s)
void GetRawDecoderInfo(art::Event &)
const double & GetCKov1Pressure() const
void SetActiveTrigger(size_t theTrigger)