BeamEvent_module.cc
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////
2 // Class: BeamEvent
3 // Plugin Type: producer (art v2_08_03)
4 // File: BeamEvent_module.cc
5 //
6 // Written by Jake Calcutt (calcuttj@msu.edu)
7 ////////////////////////////////////////////////////////////////////////
8 
9 
18 #include "fhiclcpp/ParameterSet.h"
20 #include "ifdh_art/IFBeamService/IFBeam_service.h"
21 #include "art_root_io/TFileService.h"
30 #include <bitset>
31 #include <iomanip>
32 #include <utility>
33 #include <algorithm>
34 #include <limits>
35 
36 #include "TTree.h"
37 #include "TH2F.h"
38 #include "TVectorD.h"
39 #include "TPolyMarker.h"
40 
41 namespace proto {
42  class BeamEvent;
43 }
44 
45 enum tofChan{
50 };
51 
52 typedef std::numeric_limits< double > dbl;
53 
55 public:
56  explicit BeamEvent(fhicl::ParameterSet const & p);
57 
58  // The compiler-generated destructor is fine for non-base
59  // classes without bare pointers or other resource use.
60 
61  // Plugins should not be copied or assigned.
62  BeamEvent(BeamEvent const &) = delete;
63  BeamEvent(BeamEvent &&) = delete;
64  BeamEvent & operator = (BeamEvent const &) = delete;
65  BeamEvent & operator = (BeamEvent &&) = delete;
66 
67  // Required functions.
68  //void reconfigure(fhicl::ParameterSet const & p);
69  void reset();
70  void reset_gentrigs();
71  void produce(art::Event & e) override;
72 
73  // Selected optional functions.
74  void beginJob() override;
75 
76  uint64_t joinHighLow(double,double);
77 
78  TVector3 ConvertProfCoordinates(double x, double y, double z, double zOffset);
80  bool rotated = false;
81  void RotateMonitorVector(TVector3 &vec);
82 
84 
85  void TimeIn(art::Event &, uint64_t);
86  void GetSpillInfo(art::Event &);
87  void MatchBeamToTPC();
88  void MatchS11ToGen();
89  void SetBeamEvent();
90  void MakeTrack(size_t);
91  void MomentumSpec(size_t);
92  double MomentumCosTheta(double, double, double);
93 
94  double GetPosition(std::string, int);
95  void MaskGlitches( std::vector<short> &, std::array<short,192> & );
96  TVector3 ProjectToTPC(TVector3, TVector3);
97  double GetPairedPosition(std::string, size_t);
98 
100  void parseGeneralXBPF(std::string, uint64_t, size_t);
101  void parseXBPF(uint64_t);
102 
103  void parseXTOF(uint64_t);
104  void parseXCETDB(uint64_t);
105 
106 
107  void getS11Info(uint64_t);
108 
109  std::vector<double> FetchAndReport(
110  long long, std::string,
111  std::unique_ptr<ifbeam_ns::BeamFolder>&);
112 
113 private:
114 
115  TTree * fOutTree = 0x0;
116  TH1F * fFullMomentum = 0x0;
117  TH1F * fCutMomentum = 0x0;
118 
119  TTree * fGenTrigTree = 0x0;
120  TTree * fXTOF1ATree = 0x0;
121  TTree * fXTOF1BTree = 0x0;
122  TTree * fXTOF2ATree = 0x0;
123  TTree * fXTOF2BTree = 0x0;
124 
125 
126  double fGenTrigFrac = 0;
127  double fGenTrigSec = 0;
128  double fGenTrigCoarse = 0;
129  double fXTOF1AFrac = 0;
130  double fXTOF1ACoarse = 0;
131  double fXTOF1BFrac = 0;
132  double fXTOF1BCoarse = 0;
133  double fXTOF2AFrac = 0;
134  double fXTOF2ACoarse = 0;
135  double fXTOF2BFrac = 0;
136  double fXTOF2BCoarse = 0;
137 
138  double fXTOF1ASec = 0;
139  double fXTOF1BSec = 0;
140  double fXTOF2ASec = 0;
141  double fXTOF2BSec = 0;
142  std::vector<double> diff2A;
143  std::vector<double> diff2B;
144 
145 
146  long long int eventTime = 0;
147  double SpillStart = 0.;
148  ULong_t SpillStart_alt = 0;
149  bool SpillStartValid = false;
150  bool acqStampValid = false;
151  double PrevStart=-99999.;
152  double SpillEnd = 0.;
153  double SpillOffset = 0.;
154  double ActiveTriggerTime = 0.;
155  long long RDTSTime = 0;
156  double RDTSTimeSec = 0.;
157  double PrevRDTSTimeSec=-99999.;
158  double RDTSTimeNano = 0.;
159  double RDTSSec, RDTSNano;
160 
161  //long long valid_fetch_time = 0;
162  //long long spill_valid_fetch_time = 0;
163 
164  double s11Sec = 0., s11Nano = 0.;
165  double fOutTOF= 0.;
166  double fOutP= 0.;
167  int fOutC0 = 0, fOutC1 = 0;
168 
169  std::vector<double> genTrigFracs, genTrigCoarses, genTrigSecs;
170  bool fMatched;
171 
172  int RDTSTrigger = 0;
173 
174  double acqTime = 0.;
175  double acqStampMBPL = 0.;
176 
177  int C1DB = 0;
178  int C2DB = 0;
179 
180 
181  int eventNum = 0;
182  int runNum = 0;
183  int subRunNum = 0;
184  double CKov1Pressure = 0.;
185  double CKov2Pressure = 0.;
186  int CKov1Counts = 0;
187  int CKov2Counts = 0;
188 
189  TVector3 fBMBasisX = TVector3(1.,0.,0.);
190  TVector3 fBMBasisY = TVector3(0.,1.,0.);
191  TVector3 fBMBasisZ = TVector3(0.,0.,1.);
192 
200  double fTimeWindow;
201  uint64_t fFixedTime;
202  std::vector< std::string > fDevices;
203 
208 
213  double fBeamBend;
214  double L1, L2, L3;
215  double fBProf1Shift;
216  double fBProf2Shift;
217  double fBProf3Shift;
218 
219  std::map<std::string, std::string > fDeviceTypes;
220  std::map< std::string, double > fFiberDimension;
221 
222 
223  // Names of the CERN Beam Devices
224  // Names have the form of a prefix + device
225 
226  // Prefixes for different devices classes:
227  // Beam Positions (BPF)
228  // Time of flights (TOF)
229  // and Cerenkov counters (CET)
233 
234  // Device Names for each device
235 
236  // Time of Flights
241 
243 
244  // Cerenkovs
250 
254 
257  double fNP04FrontZ;
258  double fBeamX, fBeamY, fBeamZ;
259 
265 
268  double fOffsetTAI;
269  double fS11DiffUpper;
270  double fS11DiffLower;
271 
274 
277 
280 
285 
287  bool gotTOFs;
289 
292 
295 
296  uint64_t cache_start = 0;
297  uint64_t cache_end = 0;
298 
299  std::unique_ptr<ifbeam_ns::BeamFolder> bfp;
300  std::unique_ptr<ifbeam_ns::BeamFolder> bfp_xcet;
302 
304 
305 // double L1=1.980, L2=1.69472, L3=2.11666;
306  double magnetLen = 1., magnetField = 1000.;
307  std::vector<double> current;
308 
309  //Hardware Parameters for magnetic field stuff
310  double mag_P1 = 5.82044830e-3;
311  // unused double mag_P2 = 0.;
312  double mag_P3 = -4.68880000e-6;
313  double mag_P4 = 324.573967;
314 
315 };
316 
317 // Constructor
319  : EDProducer(p),
320  fBundleName(p.get<std::string>("BundleName")),
321  fXCETBundleName(p.get<std::string>("XCETBundleName")),
322  fOutputLabel(p.get<std::string>("OutputLabel")),
323  fURLStr(p.get<std::string>("URLStr")),
324  fBFEpsilon(p.get<double>("BFEpsilon")),
325  fXCETEpsilon(p.get<double>("XCETEpsilon")),
326  fXCETFetchShift(p.get<double>("XCETFetchShift")),
327  fIFBeamDebug(p.get<int>("IFBeamDebug")),
328  fTimeWindow(p.get<double>("TimeWindow")),
329  fFixedTime(p.get<uint64_t>("FixedTime")),
330  fDevices(p.get<std::vector<std::string>>("Devices")),
331  //For Tracking/////
332  firstUpstreamName(p.get<std::string>("FirstUpstream")),
333  secondUpstreamName(p.get<std::string>("SecondUpstream")),
334  firstDownstreamName(p.get<std::string>("FirstDownstream")),
335  secondDownstreamName(p.get<std::string>("SecondDownstream")),
336  ///////////////////
337  //For Momentum Spectrometry////
338  firstBPROF1(p.get<std::string>("FirstBPROF1")),
339  secondBPROF1(p.get<std::string>("SecondBPROF1")),
340  BPROF2(p.get<std::string>("BPROF2")),
341  BPROF3(p.get<std::string>("BPROF3")),
342  fBeamBend(p.get<double>("BeamBend")),
343  L1(p.get<double>("L1")),
344  L2(p.get<double>("L2")),
345  L3(p.get<double>("L3")),
346  fBProf1Shift(p.get<double>("BProf1Shift")),
347  fBProf2Shift(p.get<double>("BProf2Shift")),
348  fBProf3Shift(p.get<double>("BProf3Shift")),
349  fXBPFPrefix(p.get<std::string>("XBPFPrefix")),
350  fXTOFPrefix(p.get<std::string>("XTOFPrefix")),
351  fXCETPrefix(p.get<std::string>("XCETPrefix")),
352  //XTOF devices
353  fTOF1(p.get< std::string >("TOF1")),
354  fTOF2(p.get< std::string >("TOF2")),
355  fTOF1A(fTOF1 + "A"),
356  fTOF1B(fTOF1 + "B"),
357  fTOF2A(fTOF2 + "A"),
358  fTOF2B(fTOF2 + "B"),
359  fTOFCalAA(p.get<double>("TOFCalAA")),
360  fTOFCalBA(p.get<double>("TOFCalBA")),
361  fTOFCalAB(p.get<double>("TOFCalAB")),
362  fTOFCalBB(p.get<double>("TOFCalBB")),
363  fCKov1(p.get< std::string >("CKov1")),
364  fCKov2(p.get< std::string >("CKov2")),
365  fXCET1(p.get< std::string >("XCET1")),
366  fXCET2(p.get< std::string >("XCET2")),
367 
368  fFillCacheUp(p.get< double >("FillCacheUp")),
369  fFillCacheDown(p.get< double >("FillCacheDown")),
370 
371 
372 
373  //New parameters to match Leigh's
374  fRotateMonitorXZ(p.get<double>("RotateMonitorXZ")),
375  fRotateMonitorYZ(p.get<double>("RotateMonitorYZ")),
376  fRotateMonitorYX(p.get<double>("RotateMonitorYX")),
377 
378  fFirstTrackingProfZ(p.get<double>("FirstTrackingProfZ")),
379  fSecondTrackingProfZ(p.get<double>("SecondTrackingProfZ")),
380  fNP04FrontZ(p.get<double>("NP04FrontZ")),
381 
382  fBeamX(p.get<double>("BeamX")),
383  fBeamY(p.get<double>("BeamY")),
384  fBeamZ(p.get<double>("BeamZ")),
385 
386  fForceNewFetch(p.get<bool>("ForceNewFetch")),
387  fXCETDebug(p.get<bool>("XCETDebug")),
388  fMatchTime(p.get<bool>("MatchTime")),
389  fForceRead(p.get<bool>("ForceRead")),
390  fForceMatchS11(p.get<bool>("ForceMatchS11")),
391 
392 
393  fTimingCalibration(p.get<double>("TimingCalibration")),
394  fCalibrationTolerance(p.get<double>("CalibrationTolerance")),
395  fOffsetTAI(p.get<double>("OffsetTAI")),
396 
397  fS11DiffUpper(p.get<double>("S11DiffUpper")),
398  fS11DiffLower(p.get<double>("S11DiffLower")),
399 
400  fRDTSToS11Upper(p.get<double>("RDTSToS11Upper")),
401  fRDTSToS11Lower(p.get<double>("RDTSToS11Lower")),
402 
403  fOffsetCTBtoRDTS(p.get<int>("OffsetCTBtoRDTS")),
404  fToleranceCTBtoRDTS(p.get<int>("ToleranceCTBtoRDTS")),
405 
406  fDownstreamToGenTrig(p.get<double>("DownstreamToGenTrig")),
407  fUpstreamToDownstream(p.get<double>("UpstreamToDownstream")),
408 
409  fPrintDebug(p.get<bool>("PrintDebug")),
410  fSaveOutTree(p.get<bool>("SaveOutTree")),
411  fDebugMomentum(p.get<bool>("DebugMomentum")),
412  fDebugTOFs(p.get<bool>("DebugTOFs"))
413 
414  {
415  // Declare products this module will provide
416  produces<std::vector<beam::ProtoDUNEBeamEvent>>();
417 
418  //magnetLen = 1.;//(m)
419  //magnetField = 1000.;//()
420  ///////////////////////////////
421 
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() );
424 
425  //Deminsion of Fibers
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());
428 
429 
430 
431 }
432 // END Constructor
433 ////////////////////////
434 
435 std::vector<double> proto::BeamEvent::FetchAndReport(long long time, std::string name, std::unique_ptr<ifbeam_ns::BeamFolder>& the_folder){
436 
437  //Note! Sometimes this won't retrieve the data from the database. We'll need
438  // catch the exception outside of this and handle it according to whichever
439  // device we're trying to retrieve from.
440 
441  std::vector<double> theResult;
442  if( fPrintDebug ){
443  MF_LOG_INFO("BeamEvent") << "Trying to grab from folder: " << name << "\n";
444  MF_LOG_INFO("BeamEvent") << "At Time: " << time << "\n";
445  }
446 
447  theResult = the_folder->GetNamedVector(time, name);
448 
449  if( fPrintDebug )
450  MF_LOG_INFO("BeamEvent") << "Successfully fetched " << time << "\n";
451 
452  return theResult;
453 }
454 
455 
456 
457 //Gets the Timing and CTB raw-decoded info.
458 //Finds the triggers, and looks for a valid trigger word
459 //(i.e. coming from beam)
460 //
461 //Returns the timestamp of the high level trigger.
463 
464  if( fPrintDebug ){
465  MF_LOG_INFO("BeamEvent") << "\n";
466  MF_LOG_INFO("BeamEvent") << "Getting Raw Decoder Info" << "\n";
467  }
468 
469  art::InputTag itag1("timingrawdecoder","daq");
470  RDTimeStampHandle = e.getHandle< std::vector<raw::RDTimeStamp> >(itag1);
471 
472  for (auto const & RDTS : *RDTimeStampHandle){
473 
474  uint64_t high = RDTS.GetTimeStamp_High();
475  uint64_t low = RDTS.GetTimeStamp_Low();
476  high = high << 32;
477  uint64_t joined = (high | low);
478 
479  RDTSTime = joined;
480  RDTSTrigger = RDTS.GetFlags();
481 
482  if( fPrintDebug ){
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";
486  MF_LOG_INFO("BeamEvent") << "Trigger: " << RDTSTrigger << "\n";
487  }
488 
489  //Separates seconds portion of the ticks
490  //From the nanoseconds
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;
494 
495  //Units are 20 nanoseconds ticks
496  RDTSTimeSec = 20.e-9 * RDTSTickSec;
497  RDTSTimeNano = 20. * RDTSTickNano;
498  RDTSSec = 20.e-9 * RDTSTickSec;
499  RDTSNano = 20. * RDTSTickNano;
500 
501  }
502 
503 }
504 
506 
507  auto const PDTStampHandle = e.getValidHandle< std::vector< dune::ProtoDUNETimeStamp > > ("timingrawdecoder:daq");
508  std::vector< dune::ProtoDUNETimeStamp > PDTStampVec(*PDTStampHandle);
509  auto PDTStamp = PDTStampVec[0];
510 
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();
515  SpillStart_alt = PDTStamp.getLastSpillStart();
516 
517 
518  if( 0ul == ~(SpillStart_alt) ) SpillStartValid = false;
519  else SpillStartValid = true;
520 
521  if( fPrintDebug ){
522  MF_LOG_INFO("BeamEvent") << PDTStamp.getLastSpillStart() << " " << SpillStart_alt << " " << SpillStartValid << "\n";
523  MF_LOG_INFO("BeamEvent") << "Version: " << ver << "\n";
524  MF_LOG_INFO("BeamEvent") << "Run: " << RunStart << "\n";
525  MF_LOG_INFO("BeamEvent") << "Spill Start: " << std::setw(20) << SpillStart << "\n";
526  MF_LOG_INFO("BeamEvent") << "Spill End: " << SpillEnd << "\n";
527  }
528 
529 
530 }
531 
533 
534  /////Now look at the acqStamp coming out of IFBeam
535  try{
536  std::vector<double> acqStamp = FetchAndReport(time, "dip/acc/NORTH/NP04/POW/MBPL022699:acqStamp[]", bfp);
537 
538  if( acqStamp[0] < 300000000.0 ){
539  if( fPrintDebug ){
540  MF_LOG_INFO("BeamEvent") << "Warning: MBPL Spill Start is low " << acqStamp[0]
541  << "\nWill need to time in with S11\n";
542  }
543  acqStampValid = false;
544  }
545  else{ acqStampValid = true; }
546 
547  acqStampMBPL = 1.e-9 * joinHighLow(acqStamp[0], acqStamp[1]);
548  if( fPrintDebug )
549  MF_LOG_INFO("BeamEvent") << "MBPL: " << acqStampMBPL << "\n";
550 
551  //Assign the calibration offset
553 
554  }
555  catch( std::exception const&){
556  acqStampValid = false;
557  MF_LOG_WARNING("BeamEvent") << "Could not get Spill time to time in\n";
558  }
559 
560 }
561 
563 
564  if( fPrintDebug )
565  MF_LOG_INFO("BeamEvent") << "Matching S11 To Gen" << "\n";
566 
567  for(size_t iT = 0; iT < beamspill->GetNT0(); ++iT){
568  double GenTrigSec = beamspill->GetT0(iT).first;
569  double GenTrigNano = beamspill->GetT0(iT).second;
570 
571  double diff = GenTrigSec - s11Sec;
572  diff += 1.e-9*(GenTrigNano - s11Nano);
573 
574  if( fS11DiffLower < diff && diff < fS11DiffUpper ){
575 
576  if( fPrintDebug ){
577  MF_LOG_INFO("BeamEvent") << "Found matching S11 and GenTrig!" << "\n";
578  MF_LOG_INFO("BeamEvent") << "diff: " << diff << "\n";
579  }
580 
581  beamspill->SetActiveTrigger( iT );
582  beamevt->SetActiveTrigger( iT );
583  beamevt->SetT0( beamspill->GetT0( iT ) );
584  return;
585  }
586  }
587 
588  if( fPrintDebug )
589  MF_LOG_INFO("BeamEvent") << "Could not find matching time " << "\n";
590 
591  beamspill->SetUnmatched();
592 }
593 
594 
596 
597  if( fPrintDebug )
598  MF_LOG_INFO("BeamEvent") << "Matching in time between Beamline and TPC!!!" << "\n";
599  for(size_t iT = 0; iT < beamspill->GetNT0(); ++iT){
600 
601  double GenTrigSec = beamspill->GetT0(iT).first;
602  double GenTrigNano = beamspill->GetT0(iT).second;
603 
604  //Separates seconds portion of the ticks
605  //From the nanoseconds
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;
609 
610  //Units are 20 nanoseconds ticks
611  double RDTSTimeSec = 20.e-9 * RDTSTickSec;
612  double RDTSTimeNano = 20. * RDTSTickNano;
613 
614  double diffSec = RDTSTimeSec - GenTrigSec - SpillOffset;
615  double diffNano = 1.e-09*(RDTSTimeNano - GenTrigNano);
616 
617  if( fPrintDebug ){
618  MF_LOG_INFO("BeamEvent") << "RDTSTimeSec - GenTrigSec " << RDTSTimeSec - GenTrigSec << "\n";
619  MF_LOG_INFO("BeamEvent") << "diff: " << diffSec << " " << diffNano << "\n";
620  }
621 
622  double diff = diffSec + diffNano;
623 
625 
626  beamspill->SetActiveTrigger( iT );
627  beamevt->SetActiveTrigger( iT );
628  beamevt->SetT0( beamspill->GetT0( iT ) );
629 
630  if( fPrintDebug ){
631  MF_LOG_INFO("BeamEvent") << "FOUND MATCHING TIME!!!" << "\n";
632  MF_LOG_INFO("BeamEvent") << "diff: " << diff << "\n";
633  MF_LOG_INFO("BeamEvent") << "Set event T0: " << beamevt->GetT0Sec() << " " << beamevt->GetT0Nano() << "\n";
634  }
635 
636  return;
637  }
638  }
639 
640  MF_LOG_INFO("BeamEvent") << "Could not find matching time " << "\n";
641  beamspill->SetUnmatched();
642 }
643 
645 
646  if( !beamspill->CheckIsMatched() ){
647  MF_LOG_INFO("BeamEvent") << "art Event is unmatched to Beam Spill " << "\n";
648  return;
649  }
650 
651  size_t activeTrigger = beamspill->GetActiveTrigger();
652  beamevt->SetActiveTrigger( activeTrigger );
653  beamevt->SetT0( beamspill->GetT0( activeTrigger ) );
654 
655  if( fPrintDebug ){
656  MF_LOG_INFO("BeamEvent") << "Setting beam event for matched event" << "\n";
657  MF_LOG_INFO("BeamEvent") << "SetActiveTrigger " << beamevt->GetActiveTrigger() << "\n";
658  MF_LOG_INFO("BeamEvent") << "Set T0 " << beamevt->GetT0Sec() << " " << beamevt->GetT0Nano() << "\n" << "\n";
659  }
660 
661 
662  for( size_t i = 0; i < fDevices.size(); ++i){
663  std::string theName = fDevices[i];
664  beamevt->SetFBMTrigger( theName, beamspill->GetFBM(theName, activeTrigger) );
665  }
666 
667  beamevt->SetTOFs( beamspill->GetMultipleTOFs( activeTrigger ) );
668  beamevt->SetTOFChans( beamspill->GetMultipleTOFChans( activeTrigger ) );
669  beamevt->SetUpstreamTriggers( beamspill->GetUpstreamTriggers( activeTrigger ) );
670  beamevt->SetDownstreamTriggers( beamspill->GetDownstreamTriggers( activeTrigger ) );
673  beamevt->DecodeTOF();
674  beamevt->SetMagnetCurrent( beamspill->GetMagnetCurrent() );
675  beamevt->SetCKov0( beamspill->GetCKov0( activeTrigger ) );
676  beamevt->SetCKov1( beamspill->GetCKov1( activeTrigger ) );
677 
678 
679  if( fPrintDebug ){
680  MF_LOG_INFO("BeamEvent") << "beamevt has TOF " << beamevt->GetTOF()
681  << " and TOFChan " << beamevt->GetTOFChan() << "\n";
682 
683 
684  MF_LOG_INFO("BeamEvent") << "beamevt has Magnet Current " << beamevt->GetMagnetCurrent() << "\n";
685  MF_LOG_INFO("BeamEvent") << "beamevt CKov0: " << beamevt->GetCKov0Status() << " " << beamevt->GetCKov0Pressure() << "\n";
686  MF_LOG_INFO("BeamEvent") << "beamevt CKov1: " << beamevt->GetCKov1Status() << " " << beamevt->GetCKov1Pressure() << "\n";
687  MF_LOG_INFO("BeamEvent") << "Finished adding info to beamevt " << "\n";
688 
689  }
690 
691  C1DB = beamspill->GetCKov0Status( activeTrigger );
692  C2DB = beamspill->GetCKov1Status( activeTrigger );
693 }
694 
695 
697  genTrigFracs.clear();
698  genTrigCoarses.clear();
699  genTrigSecs.clear();
700 }
701 
703  acqTime = 0;
704  acqStampMBPL = 0;
705  RDTSTrigger = -1;
706  C1DB = -1;
707  C2DB = -1;
708  SpillStart = -1;
709  SpillEnd = -1;
710  SpillOffset = -1;
711 
712  s11Nano = -1.;
713  s11Sec = -1.;
714 
715  ActiveTriggerTime = -1;
716  RDTSTime = 0;
717  RDTSSec = -999.;
718  RDTSNano = -999.;
719 
720  fOutP = -1.;
721  fOutTOF = -1.;
722  fOutC0 = -1;
723  fOutC1 = -1;
724  fMatched = false;
725 }
726 
727 
728 ////////////////////////
729 // Producer Method (reads in the event and derives values)
731 
732  reset();
733 
734  eventNum = e.event();
735  runNum = e.run();
736  subRunNum = e.subRun();
737 
738  if( e.time().timeHigh() == 0 ) eventTime = e.time().timeLow();
739  else eventTime = e.time().timeHigh();
740 
741 
742 
743  // Create a new beam event (note the "new" here)
746 
747  //Get the information from the timing system
748  //This gets RDTSTrigger and RDTSTime
750 
751  //Get Spill Information
752  //This stores Spill Start, Spill End,
753  //And Prev Spill Start
754  GetSpillInfo(e);
755 
756 
757  //Check if we have a valid beam trigger
758  //If not, just place an empty beamevt
759  //and move on
760  //
761  //Also check if we've gotten good spill info
762  //
763  //Or if we're forcing to read out the Beamline Info
764  if( ( (RDTSTrigger == 12) ) || fForceRead ){
765 
766  //Start getting beam event info
767  uint64_t fetch_time = uint64_t( RDTSTime * 2e-8 );
768  uint64_t fetch_time_down = uint64_t( RDTSTime * 2e-8 );
769 
770  if( fPrintDebug )
771  MF_LOG_INFO("BeamEvent") << "RDTSTime: " << uint64_t( RDTSTime * 2e-8 ) << "\n";
772 
773  //Check if we are still using the same spill information.
774  //
775  //If it's a new spill: Get new info from the database
776  //Each 'parse' command fetches new data from the database
777  //
778  //If it's the same spill then just pass the old BeamEvent
779  //Object. Its 'active trigger' info will be updated below
780  //
781  //Also: Check if the SpillStart is invalid. If the RDTSTime is
782  //Greater than 5, it's in a new spill, so fetch again
783  //
784  //Can be overridden with a flag from the fcl
785  if( ( ( PrevStart != SpillStart) || ( !SpillStartValid && ( abs(RDTSTimeSec - PrevRDTSTimeSec) > 5 ) ) )
786  || fForceNewFetch){
787  reset_gentrigs();
788 
789  if( fPrintDebug )
790  MF_LOG_INFO("BeamEvent") << "New spill or forced new fetch. Getting new beamspill info" << "\n";
791 
792  //Testing: printing out cache start and end
793  cache_start = bfp->GetCacheStartTime();
794  cache_end = bfp->GetCacheEndTime();
795 
796  if( fPrintDebug ){
797  MF_LOG_INFO("BeamEvent") << "cache_start: " << cache_start << "\n";
798  MF_LOG_INFO("BeamEvent") << "cache_end: " << cache_end << "\n";
799  MF_LOG_INFO("BeamEvent") << "fetch_time: " << fetch_time << "\n";
800  }
801 
802  cache_start = bfp_xcet->GetCacheStartTime();
803  cache_end = bfp_xcet->GetCacheEndTime();
804 
805  if( fPrintDebug ){
806  MF_LOG_INFO("BeamEvent") << "xcet cache_start: " << cache_start << "\n";
807  MF_LOG_INFO("BeamEvent") << "xcet cache_end: " << cache_end << "\n";
808  MF_LOG_INFO("BeamEvent") << "xcet fetch_time: " << fetch_time << "\n";
809  }
810 
811  //Not the first event
812  if(cache_start > 0 && cache_end > 0){
813 
814  //So try filling the cache first with the 'possible' end of spill time
815  //then the lower spill time.
816  //
817  //This is done so that the cache essentially reshuffles where it starts and ends
818  //
819  //All the checking is done internal to the FillCache method
820  //
821  //Note: I'm using a loose definition of the start and end of spills
822  // It's really just the maximum and minimum possible vales of those times
823  // since it's not possible to know for certain in any given event.
824  // (The info does exist in the raw decoder info, but it's not always
825  // present, so I'm just opting for this)
826  try{
827  bfp->FillCache( fetch_time + fFillCacheUp );
828 
829  if( fPrintDebug ){
830  cache_start = bfp->GetCacheStartTime();
831  cache_end = bfp->GetCacheEndTime();
832  MF_LOG_INFO("BeamEvent") << "interim cache_start: " << cache_start << "\n";
833  MF_LOG_INFO("BeamEvent") << "interim cache_end: " << cache_end << "\n";
834  }
835 
836  bfp->FillCache( fetch_time - fFillCacheDown );
837  if( fPrintDebug ){
838  cache_start = bfp->GetCacheStartTime();
839  cache_end = bfp->GetCacheEndTime();
840  MF_LOG_INFO("BeamEvent") << "new cache_start: " << cache_start << "\n";
841  MF_LOG_INFO("BeamEvent") << "new cache_end: " << cache_end << "\n";
842  }
843 
844  bfp_xcet->FillCache( fetch_time + fFillCacheUp );
845  if( fPrintDebug ){
846  cache_start = bfp_xcet->GetCacheStartTime();
847  cache_end = bfp_xcet->GetCacheEndTime();
848  MF_LOG_INFO("BeamEvent") << "interim xcet cache_start: " << cache_start << "\n";
849  MF_LOG_INFO("BeamEvent") << "interim xcet cache_end: " << cache_end << "\n";
850  }
851 
852  bfp_xcet->FillCache( fetch_time - fFillCacheDown );
853  if( fPrintDebug ){
854  cache_start = bfp_xcet->GetCacheStartTime();
855  cache_end = bfp_xcet->GetCacheEndTime();
856  MF_LOG_INFO("BeamEvent") << "new xcet cache_start: " << cache_start << "\n";
857  MF_LOG_INFO("BeamEvent") << "new xcet cache_end: " << cache_end << "\n";
858  }
859  }
860  catch( std::exception const& e){
861  MF_LOG_WARNING("BeamEvent") << "Could not fill cache\n";
862  MF_LOG_ERROR("BeamEvent") << e.what() << "\n";
863  }
864  }
865  else{
866  //First event, let's get the start of spill info
867 
868  if( fPrintDebug )
869  MF_LOG_INFO("BeamEvent") << "First Event: Priming cache\n";
870 
871  try{
872  bfp->FillCache( fetch_time - fFillCacheDown );
873  }
874  catch( std::exception const& e){
875  MF_LOG_WARNING("BeamEvent") << "Could not fill cache\n";
876  MF_LOG_ERROR("BeamEvent") << e.what() << "\n";
877  }
878  try{
879  bfp_xcet->FillCache( fetch_time - fFillCacheDown );
880  }
881  catch( std::exception const& e){
882  MF_LOG_WARNING("BeamEvent") << "Could not fill xcet cache\n";
883  MF_LOG_ERROR("BeamEvent") << e.what() << "\n";
884  }
885 
886  if( fPrintDebug ){
887  cache_start = bfp->GetCacheStartTime();
888  cache_end = bfp->GetCacheEndTime();
889  MF_LOG_INFO("BeamEvent") << "new cache_start: " << cache_start << "\n";
890  MF_LOG_INFO("BeamEvent") << "new cache_end: " << cache_end << "\n";
891 
892  cache_start = bfp_xcet->GetCacheStartTime();
893  cache_end = bfp_xcet->GetCacheEndTime();
894  MF_LOG_INFO("BeamEvent") << "new xcet cache_start: " << cache_start << "\n";
895  MF_LOG_INFO("BeamEvent") << "new xcet cache_end: " << cache_end << "\n";
896  }
897 
898  }
899 
900  // Parse the Time of Flight Counter data for the list
901  // of times that we are using
902  parseXTOF(fetch_time);
903 
904 
905  // Parse the Beam postion counter information for the list
906  // of time that we are using
908 
909  if( gotGeneralTrigger ){
910  parseXBPF(fetch_time);
911  /*
912  for(size_t d = 0; d < fDevices.size(); ++d){
913  std::string name = fDevices[d];
914  parseGeneralXBPF(name, time, d);
915  }
916  */
917  }
918 
919  parseXCETDB(fetch_time);
920 
921  try{
922  current = FetchAndReport(fetch_time_down, "dip/acc/NORTH/NP04/POW/MBPL022699:current", bfp);
923  gotCurrent = true;
924  beamspill->SetMagnetCurrent(current[0]);
925  }
926  catch( std::exception const&){
927  MF_LOG_WARNING("BeamEvent") << "Could not get magnet current\n";
928  gotCurrent = false;
929  }
930 
931  //Set PrevStart to SpillStart here
932  //so that we don't skip in the case
933  //the first event in the spill did not
934  //have a good beamline trigger
937 
938  }
939  else{
940  if( fPrintDebug )
941  MF_LOG_INFO("BeamEvent") << "Same spill. Reusing beamspill info" << "\n";
942 
944  }
945 
946  if( fPrintDebug ){
947  MF_LOG_INFO("BeamEvent") << "NGoodParticles: " << beamspill->GetNT0() << "\n";
948  MF_LOG_INFO("BeamEvent") << "NTOF0: " << beamspill->GetNTOF0Triggers() << "\n";
949  MF_LOG_INFO("BeamEvent") << "NTOF1: " << beamspill->GetNTOF1Triggers() << "\n";
950  MF_LOG_INFO("BeamEvent") << "acqTime: " << acqTime << "\n";
951  MF_LOG_INFO("BeamEvent") << "NXBPF: " << beamspill->GetNFBMTriggers(fDevices[0]) << "\n";
952  }
953 
954  if( fMatchTime ){
955 
956  //Now do the matching in Time:
957  //Get the conversion from the ProtoDUNE Timing system
958  //To the one in the SPS.
959  //
960 
962  TimeIn(e, fetch_time_down);
963 
964  if( fPrintDebug )
965  MF_LOG_INFO("BeamEvent") << "SpillOffset " << SpillOffset << "\n";
966 
967  //If not successfully timed in,
968  //Oh well. It won't cause a crash
969  //But it won't be matched.
970  //
971  //Additionally, check if the MBPL timestamp is valid,
972  //There are some instances in the database of it being abnormally low
973  //If so, then just do the S11 matching
974  if( acqStampValid ){
975  MatchBeamToTPC();
976  }
977  else{
978  try{
979  getS11Info(fetch_time);
980  }
981  catch( std::exception const& ){
982  MF_LOG_WARNING("BeamEvent") << "Could not get S11 Info\n";
983  }
984 
985  //Again, it won't crash, but it won't match
986  MatchS11ToGen();
987  }
988  }
989  else{
990  try{
991  getS11Info(fetch_time);
992  }
993  catch( std::exception const& ){
994  MF_LOG_WARNING("BeamEvent") << "Could not get S11 Info\n";
995  }
996 
997  //Again, it won't crash, but it won't match
998  MatchS11ToGen();
999  }
1000 
1001 
1002  if( beamspill->CheckIsMatched() ){
1003  std::pair<double,double> theTime = beamspill->GetT0(beamspill->GetActiveTrigger());
1004  ActiveTriggerTime = theTime.first + theTime.second*1.e-9;
1005 
1006  if( fPrintDebug )
1007  MF_LOG_INFO("BeamEvent") << "Trigger: " << beamspill->GetActiveTrigger() << " " << ActiveTriggerTime << "\n";
1008 
1009  //Pass the information to the beamevent
1010  SetBeamEvent();
1011 
1012  MakeTrack( beamspill->GetActiveTrigger() );
1013 
1014  if( fPrintDebug )
1015  MF_LOG_INFO("BeamEvent") << "Added " << beamevt->GetNBeamTracks() << " tracks to the beam spill" << "\n";
1016 
1017  //Momentum
1018  if( gotCurrent ){
1019  MomentumSpec( beamspill->GetActiveTrigger() );
1020  }
1021 
1022  if( fPrintDebug )
1023  MF_LOG_INFO("BeamEvent") << "Got NRecoBeamMomenta: " << beamevt->GetNRecoBeamMomenta() << "\n";
1024  }
1025  }
1026 
1027 
1028  //Pass beamspill to the next event;
1029  //Erase the Track and Reco Momentum info
1034 
1035  }
1036  //Start of a new spill, but the first event was
1037  //Not a beam trigger. In this case, it would not
1038  //have been filled with info in the block above
1039  //So let's make it empty so we aren't putting
1040  //old spill info in the new event
1041 
1042  //Or. If this was not a beam trigger, and the RDTSTime
1043  //Comes from a new spill while the SpillStart was invalid, pass it.
1044  else if( ( ( PrevStart != SpillStart ) || ( !SpillStartValid && ( abs(RDTSTimeSec - PrevRDTSTimeSec) > 5 ) ) )
1045  && RDTSTrigger != 12 ){
1047  }
1048 
1049  //Can Remove BITrigger,CTBTimestamp
1050  beamevt->SetBITrigger(-1);
1054  beamevt->SetCTBTimestamp( -1. );
1056 
1057  std::unique_ptr<std::vector<beam::ProtoDUNEBeamEvent> > beamData(new std::vector<beam::ProtoDUNEBeamEvent>);
1058  beamData->push_back(beam::ProtoDUNEBeamEvent(*beamevt));
1059  e.put(std::move(beamData));
1060 
1061  // Write out the to tree
1062  if( fSaveOutTree ){
1063 
1064  std::cout << "Timing trigger: " << beamevt->GetTimingTrigger() << std::endl;
1065  std::cout << "Matched: " << beamevt->CheckIsMatched() << std::endl;
1066 
1067  if( beamevt->GetTOFChan() > -1 ){
1068  fOutTOF = beamevt->GetTOF();
1069  }
1070  else fOutTOF = -1.;
1071  if( beamevt->GetNRecoBeamMomenta() == 1 ){
1073  }
1074  else fOutP = -1.;
1077  std::cout << "CKovs: " << beamevt->GetCKov0Status() << " " << beamevt->GetCKov1Status() << std::endl;
1078  std::cout << "TOF, P: " << fOutTOF << " " << fOutP << std::endl;
1079 
1080 
1081 
1083  fOutTree->Fill();
1084  }
1085 
1086  delete beamevt;
1087  delete beamspill;
1088 
1089 
1090 }
1091 // END BeamEvent::produce
1092 ////////////////////////
1093 
1095  // Places a dummy trigger vector for each device
1096 
1097  // Make a vector the names of each of the devices being readout
1098  std::vector<std::string> monitors;
1099  size_t nDev = 0;
1100 
1101  for(size_t id = 0; id < fDevices.size(); ++id){
1102  std::string name = fDevices[id];
1103  // Put the current device name on the list
1104  monitors.push_back(name);
1105  nDev++;
1106  }
1107 
1108  beamspill->InitFBMs(monitors);
1109 }
1110 // END BeamEvent::InitXBPFInfo
1111 ////////////////////////
1112 
1114  MF_LOG_INFO("BeamEvent") << "Getting S11 Info " << "\n";
1115 
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);
1120 
1121  int s11Count = (int)timestampCountS11[0];
1122 
1123  if( fPrintDebug )
1124  MF_LOG_INFO("BeamEvent") << "Found " << s11Count << " returned signals" << "\n";
1125 
1126  //Separates seconds portion of the ticks
1127  //From the nanoseconds
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;
1131 
1132  //Units are 20 nanoseconds ticks
1133  double RDTSTimeSec = 20.e-9 * RDTSTickSec;
1134  double RDTSTimeNano = 20. * RDTSTickNano;
1135 
1136 
1137  for(int i = 0; i < s11Count; ++i){
1138  double nano = 8.*coarseS11[i] + fracS11[i]/512.;
1139 
1140  double diffSec = secondsS11[2*i + 1] - RDTSTimeSec - fOffsetTAI;
1141  double diffNano = nano - RDTSTimeNano;
1142 
1143  if( fPrintDebug ){
1144  MF_LOG_INFO("BeamEvent") << i << " diffSec " << diffSec << "\n";
1145  MF_LOG_INFO("BeamEvent") << i << " diffNano " << diffNano << "\n";
1146  }
1147 
1148  double diff = diffSec + 1.e-9*diffNano;
1149  if( fRDTSToS11Lower < diff && diff < fRDTSToS11Upper ){
1150 
1151  if( fPrintDebug )
1152  MF_LOG_INFO("BeamEvent") << "FOUND Match between S11 and RDTS" << "\n";
1153 
1154  s11Nano = nano;
1155  s11Sec = secondsS11[2*i + 1] - fOffsetTAI;
1156  return;
1157  }
1158  }
1159 
1160  MF_LOG_WARNING("BeamEvent") << "Could not match RDTS to S11\n";
1161 
1162 }
1163 
1164 
1166 
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;
1172 
1173 
1174  try{
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);
1180 
1181  if( fPrintDebug )
1182  MF_LOG_INFO("BeamEvent") << "timestampCounts: " << timestampCountGeneralTrigger[0] << "\n";
1183 
1184  gotGeneralTrigger = true;
1185 
1186  uint64_t low = (uint64_t)acqStampGeneralTrigger[1];
1187  uint64_t high = (uint64_t)acqStampGeneralTrigger[0];
1188  high = high << 32;
1189  acqTime = ( high | low ) / 1000000000.;
1190 
1191  }
1192  catch(std::exception const&){
1193  MF_LOG_WARNING("BeamEvent") << "Could not get GeneralTrigger information!!" << "\n";
1194  gotGeneralTrigger = false;
1195  return;
1196  }
1197 
1198  std::vector<double> coarseTOF1A;
1199  std::vector<double> fracTOF1A;
1200  std::vector<double> secondsTOF1A;
1201 
1202  std::vector<double> coarseTOF1B;
1203  std::vector<double> fracTOF1B;
1204  std::vector<double> secondsTOF1B;
1205 
1206  std::vector<double> coarseTOF2A;
1207  std::vector<double> fracTOF2A;
1208  std::vector<double> secondsTOF2A;
1209 
1210  std::vector<double> coarseTOF2B;
1211  std::vector<double> fracTOF2B;
1212  std::vector<double> secondsTOF2B;
1213 
1214  try{
1215  coarseTOF1A = FetchAndReport(time, fXTOFPrefix + fTOF1A + ":coarse[]", bfp);
1216  fracTOF1A = FetchAndReport(time, fXTOFPrefix + fTOF1A + ":frac[]", bfp);
1217  secondsTOF1A = FetchAndReport(time, fXTOFPrefix + fTOF1A + ":seconds[]", bfp);
1218 
1219  coarseTOF1B = FetchAndReport(time, fXTOFPrefix + fTOF1B + ":coarse[]", bfp);
1220  fracTOF1B = FetchAndReport(time, fXTOFPrefix + fTOF1B + ":frac[]", bfp);
1221  secondsTOF1B = FetchAndReport(time, fXTOFPrefix + fTOF1B + ":seconds[]", bfp);
1222 
1223  coarseTOF2A = FetchAndReport(time, fXTOFPrefix + fTOF2A + ":coarse[]", bfp);
1224  fracTOF2A = FetchAndReport(time, fXTOFPrefix + fTOF2A + ":frac[]", bfp);
1225  secondsTOF2A = FetchAndReport(time, fXTOFPrefix + fTOF2A + ":seconds[]", bfp);
1226 
1227  coarseTOF2B = FetchAndReport(time, fXTOFPrefix + fTOF2B + ":coarse[]", bfp);
1228  fracTOF2B = FetchAndReport(time, fXTOFPrefix + fTOF2B + ":frac[]", bfp);
1229  secondsTOF2B = FetchAndReport(time, fXTOFPrefix + fTOF2B + ":seconds[]", bfp);
1230 
1231  gotTOFs = true;
1232  }
1233  catch(std::exception const&){
1234  MF_LOG_WARNING("BeamEvent") << "Could not get TOF information!!" << "\n";
1235  gotTOFs = false;
1236  }
1237 
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;
1243 
1244 
1245  for(size_t i = 0; i < timestampCountGeneralTrigger[0]; ++i){
1246 
1247  if( fPrintDebug ){
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";
1250  }
1251 
1252  //2*i + 1 because the format is weird
1253  fGenTrigSec = secondsGeneralTrigger[2*i + 1];
1254  fGenTrigCoarse = coarseGeneralTrigger[i];
1255  fGenTrigFrac = fracGeneralTrigger[i];
1256 
1257  genTrigSecs.push_back(fGenTrigSec);
1258  genTrigFracs.push_back(fGenTrigFrac);
1259  genTrigCoarses.push_back(fGenTrigCoarse);
1260 
1261  unorderedGenTrigTime.push_back( std::make_pair(fGenTrigSec, (fGenTrigCoarse*8. + fGenTrigFrac/512.)) );
1262 
1263  if (fGenTrigCoarse == 0.0 && fGenTrigFrac == 0.0 && fGenTrigSec == 0.0) break;
1264  if(fDebugTOFs)fGenTrigTree->Fill();
1265  }
1266 
1267  if( gotTOFs ){
1268 
1269  for(size_t i = 0; i < timestampCountGeneralTrigger[0]; ++i){
1270 
1271  if( fPrintDebug )
1272  MF_LOG_INFO("BeamEvent") << "TOF1A " << i << " " << secondsTOF1A[2*i+1] << " " << 8.*coarseTOF1A[i] + fracTOF1A[i]/512. << "\n";
1273 
1274  fXTOF1ACoarse = coarseTOF1A[i];
1275  fXTOF1AFrac = fracTOF1A[i];
1276  fXTOF1ASec = secondsTOF1A[2*i + 1];
1277 
1278  if(fXTOF1ASec < secondsTOF1A[1]) break;
1279 
1280  if (fXTOF1ACoarse == 0.0 && fXTOF1AFrac == 0.0 && fXTOF1ASec == 0.0) break;
1281  unorderedTOF1ATime.push_back(std::make_pair(fXTOF1ASec, (fXTOF1ACoarse*8. + fXTOF1AFrac/512.)) );
1282 
1283  if(fDebugTOFs)fXTOF1ATree->Fill();
1284  }
1285 
1286  for(size_t i = 0; i < timestampCountGeneralTrigger[0]; ++i){
1287 
1288  if( fPrintDebug )
1289  MF_LOG_INFO("BeamEvent") << "TOF1B " << i << " " << secondsTOF1B[2*i+1] << " " << 8.*coarseTOF1B[i] + fracTOF1B[i]/512. << "\n";
1290 
1291  fXTOF1BCoarse = coarseTOF1B[i];
1292  fXTOF1BFrac = fracTOF1B[i];
1293  fXTOF1BSec = secondsTOF1B[2*i + 1];
1294 
1295  if(fXTOF1BSec < secondsTOF1B[1]) break;
1296 
1297  if (fXTOF1BCoarse == 0.0 && fXTOF1BFrac == 0.0 && fXTOF1BSec == 0.0) break;
1298  unorderedTOF1BTime.push_back(std::make_pair(fXTOF1BSec, (fXTOF1BCoarse*8. + fXTOF1BFrac/512.)) );
1299  if(fDebugTOFs)fXTOF1BTree->Fill();
1300  }
1301 
1302  for(size_t i = 0; i < timestampCountGeneralTrigger[0]; ++i){
1303 
1304  if( fPrintDebug )
1305  MF_LOG_INFO("BeamEvent") << "TOF2A " << i << " " << secondsTOF2A[2*i+1] << " " << 8.*coarseTOF2A[i] + fracTOF2A[i]/512. << "\n";
1306 
1307  fXTOF2ACoarse = coarseTOF2A[i];
1308  fXTOF2AFrac = fracTOF2A[i];
1309  fXTOF2ASec = secondsTOF2A[2*i + 1];
1310 
1311  if(fXTOF2ASec < secondsTOF2A[1]) break;
1312 
1313  if (fXTOF2ACoarse == 0.0 && fXTOF2AFrac == 0.0 && fXTOF2ASec == 0.0) break;
1314  unorderedTOF2ATime.push_back(std::make_pair(fXTOF2ASec, (fXTOF2ACoarse*8. + fXTOF2AFrac/512.)) );
1315 
1316  if(diff2A.size()) diff2A.clear();
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) );
1319  }
1320 
1321  if(fDebugTOFs)fXTOF2ATree->Fill();
1322 
1323  }
1324 
1325  for(size_t i = 0; i < timestampCountGeneralTrigger[0]; ++i){
1326 
1327  if( fPrintDebug )
1328  MF_LOG_INFO("BeamEvent") << "TOF2B " << i << " " << secondsTOF2B[2*i+1] << " " << 8.*coarseTOF2B[i] + fracTOF2B[i]/512. << "\n";
1329 
1330  fXTOF2BCoarse = coarseTOF2B[i];
1331  fXTOF2BFrac = fracTOF2B[i];
1332  fXTOF2BSec = secondsTOF2B[2*i + 1];
1333 
1334  if(fXTOF2BSec < secondsTOF2B[1]) break;
1335  if (fXTOF2BCoarse == 0.0 && fXTOF2BFrac == 0.0 && fXTOF2BSec == 0.0) break;
1336 
1337  unorderedTOF2BTime.push_back(std::make_pair(fXTOF2BSec, (fXTOF2BCoarse*8. + fXTOF2BFrac/512.) ));
1338 
1339  if(diff2B.size()) diff2B.clear();
1340 
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) );
1343  }
1344 
1345  if(fDebugTOFs)fXTOF2BTree->Fill();
1346  }
1347  }
1348 
1349  if( fPrintDebug )
1350  MF_LOG_INFO("BeamEvent") << "NGenTrigs: " << timestampCountGeneralTrigger[0] << " NTOF2s: " << unorderedTOF2ATime.size() + unorderedTOF2BTime.size() << "\n";
1351 
1352  for(size_t iT = 0; iT < unorderedGenTrigTime.size(); ++iT){
1353 
1354  bool found_TOF = false;
1355 
1356  double the_gen_sec = unorderedGenTrigTime[iT].first;
1357  double the_gen_ns = unorderedGenTrigTime[iT].second;
1358 
1359  //1A2A = 0; 1B2A = 1, 1A2B = 2, 1B2B = 3
1360  //Add 1 for 1B, add 2 for 2B
1361  int channel = 0;
1362 
1363  std::vector< double > possibleTOF;
1364  std::vector< int > possibleTOFChan;
1365  std::vector< size_t > UpstreamTriggers;
1366  std::vector< size_t > DownstreamTriggers;
1367 
1368  if( gotTOFs ){
1369 
1370  if( fPrintDebug )
1371  MF_LOG_INFO("BeamEvent") << "Gen: " << the_gen_sec << " " << the_gen_ns << "\n";
1372 
1373  //First check 2A
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;
1378 
1379  if( delta_2A < 0. ) break;
1380  else if( delta_2A > fDownstreamToGenTrig ) continue;
1381 
1382  //if here, 0. < delta < 50ns
1383  if( fPrintDebug )
1384  MF_LOG_INFO("BeamEvent") << "Found match 2A to Gen" << "\n";
1385 
1386  //So check 1A and 1B
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;
1391 
1392  if( delta < 0. ) break;
1393  else if( delta > 0. && delta < fUpstreamToDownstream){
1394  if( fPrintDebug )
1395  MF_LOG_INFO("BeamEvent") << "Found match 1A to 2A " << delta << "\n";
1396 
1397  found_TOF = true;
1398  channel = k1A2A;
1399 
1400  possibleTOF.push_back( delta );
1401  possibleTOFChan.push_back( channel );
1402 
1403  UpstreamTriggers.push_back( ip1A );
1404  DownstreamTriggers.push_back( ip2A );
1405  }
1406  else continue;
1407  }
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;
1412 
1413  if( delta < 0. ) break;
1414  else if( delta > 0. && delta < fUpstreamToDownstream){
1415 
1416  if( fPrintDebug )
1417  MF_LOG_INFO("BeamEvent") << "Found match 1B to 2A " << delta << "\n";
1418 
1419  found_TOF = true;
1420  channel = k1B2A;
1421 
1422  possibleTOF.push_back( delta );
1423  possibleTOFChan.push_back( channel );
1424 
1425  UpstreamTriggers.push_back( ip1B );
1426  DownstreamTriggers.push_back( ip2A );
1427  }
1428  else continue;
1429  }
1430  }
1431 
1432  //Then check 2B
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;
1437 
1438 
1439  if( delta_2B < 0. ) break;
1440  else if( delta_2B > fDownstreamToGenTrig ) continue;
1441 
1442  //if here, 0. < delta < 50ns
1443  if( fPrintDebug )
1444  MF_LOG_INFO("BeamEvent") << "Found match 2B to Gen" << "\n";
1445 
1446  //So check 1A and 1B
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;
1451 
1452 
1453  if( delta < 0. ) break;
1454  else if( delta > 0. && delta < fUpstreamToDownstream){
1455 
1456  if( fPrintDebug )
1457  MF_LOG_INFO("BeamEvent") << "Found match 1A to 2B " << delta << "\n";
1458 
1459  found_TOF = true;
1460  channel = k1A2B;
1461 
1462  possibleTOF.push_back( delta );
1463  possibleTOFChan.push_back( channel );
1464 
1465  UpstreamTriggers.push_back( ip1A );
1466  DownstreamTriggers.push_back( ip2B );
1467  }
1468  else continue;
1469  }
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;
1474 
1475 
1476  if( delta < 0. ) break;
1477  else if( delta > 0. && delta < fUpstreamToDownstream){
1478 
1479  if( fPrintDebug )
1480  MF_LOG_INFO("BeamEvent") << "Found match 1B to 2B " << delta << "\n";
1481 
1482  found_TOF = true;
1483  channel = k1B2B;
1484 
1485  possibleTOF.push_back( delta );
1486  possibleTOFChan.push_back( channel );
1487 
1488  UpstreamTriggers.push_back( ip1B );
1489  DownstreamTriggers.push_back( ip2B );
1490  }
1491  else continue;
1492  }
1493  }
1494 
1495  if( fPrintDebug )
1496  MF_LOG_INFO("BeamEvent") << "Found " << possibleTOF.size() << " matched TOFs" << "\n";
1497  }
1498 
1499  if( !found_TOF ){
1500 
1501  if( fPrintDebug )
1502  MF_LOG_INFO("BeamEvent") << "No matching TOFs found. Placing dummy\n";
1503 
1504  possibleTOF.push_back( 0. );
1505  possibleTOFChan.push_back( -1 );
1506  UpstreamTriggers.push_back(0);
1507  DownstreamTriggers.push_back(0);
1508  }
1509 
1510  beamspill->AddT0(std::make_pair(the_gen_sec - fOffsetTAI, the_gen_ns));
1511  beamspill->AddMultipleTOFs( possibleTOF );
1512  beamspill->AddMultipleTOFChans( possibleTOFChan );
1513  beamspill->AddUpstreamTriggers( UpstreamTriggers );
1514  beamspill->AddDownstreamTriggers( DownstreamTriggers );
1515 
1516  }
1517 
1518 }
1519 // END BeamEvent::parseXTOF
1520 ////////////////////////
1521 
1523 
1524  if(fCKov1 != ""){
1525 
1526  try{
1527  std::vector< double > pressureCKov1 = FetchAndReport(time, fXCETPrefix + fCKov1 + ":pressure", bfp);
1528  std::vector< double > countsCKov1 = FetchAndReport(time, fXCETPrefix + fCKov1 + ":counts", bfp);
1529 
1530  CKov1Pressure = pressureCKov1[0];
1531  CKov1Counts = countsCKov1[0];
1532  }
1533  catch( std::exception const&){
1534  MF_LOG_WARNING("BeamEvent") << "Could not get Cerenkov 1 Pressure\n";
1535  CKov1Pressure = 0.;
1536  CKov1Counts = 0.;
1537  }
1538 
1539  }
1540 
1541  if(fCKov2 != ""){
1542  try{
1543  std::vector< double > pressureCKov2 = FetchAndReport(time, fXCETPrefix + fCKov2 + ":pressure", bfp);
1544  std::vector< double > countsCKov2 = FetchAndReport(time, fXCETPrefix + fCKov2 + ":counts", bfp);
1545  CKov2Pressure = pressureCKov2[0];
1546  CKov2Counts = countsCKov2[0];
1547  }
1548  catch( std::exception const&){
1549  MF_LOG_WARNING("BeamEvent") << "Could not get Cerenkov 2 Pressure\n";
1550  CKov2Pressure = 0.;
1551  CKov2Counts = 0.;
1552  }
1553  }
1554 
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;
1559 
1560  bool fetched_XCET1=false;
1561  bool fetched_XCET2=false;
1562  if( fXCET1 != "" ){
1563  try{
1564  XCET1_seconds = FetchAndReport( time + fXCETFetchShift, fXCET1 + ":SECONDS" , bfp_xcet);
1565  XCET1_frac = FetchAndReport( time + fXCETFetchShift, fXCET1 + ":FRAC" , bfp_xcet);
1566  XCET1_coarse = FetchAndReport( time + fXCETFetchShift, fXCET1 + ":COARSE" , bfp_xcet);
1567  fetched_XCET1 = true;
1568 
1569  if( fXCETDebug ){
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;
1572  }
1573  }
1574 
1575  }
1576  catch( const std::exception &e ){
1577  MF_LOG_WARNING("BeamEvent") << "Could not get XCET1 info\n";
1578  fetched_XCET1 = false;
1579  }
1580  }
1581  if( fXCET2 != "" ){
1582  try{
1583  XCET2_seconds = FetchAndReport( time + fXCETFetchShift, fXCET2 + ":SECONDS" , bfp_xcet);
1584  XCET2_frac = FetchAndReport( time + fXCETFetchShift, fXCET2 + ":FRAC" , bfp_xcet);
1585  XCET2_coarse = FetchAndReport( time + fXCETFetchShift, fXCET2 + ":COARSE" , bfp_xcet);
1586  fetched_XCET2 = true;
1587 
1588  if( fXCETDebug ){
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;
1591  }
1592  }
1593  }
1594  catch( const std::exception &e ){
1595  MF_LOG_WARNING("BeamEvent") << "Could not get XCET2 info\n";
1596 
1597  fetched_XCET2 = false;
1598  }
1599  }
1600 
1601 
1602 
1603  //Go through the general triggers and try to match. If one can't be found, then just add a 0
1604  for( size_t i = 0; i < beamspill->GetNT0(); ++i ){
1605 
1606  if( fXCETDebug ) std::cout << "GenTrig: " << i << " " << beamspill->GetT0Sec(i) << " " << beamspill->GetT0Nano(i) << std::endl;
1607 
1608  beam::CKov status_1;
1609  status_1.pressure = CKov1Pressure;
1610  status_1.timeStamp = -1.;
1611 
1612  if( !fetched_XCET1 ) status_1.trigger = -1;
1613  else{
1614 
1615  status_1.trigger = 0;
1616 
1617  for( size_t ic1 = 0; ic1 < XCET1_seconds.size(); ++ic1 ){
1618  double delta = 1.e9 * ( beamspill->GetT0Sec(i) - (XCET1_seconds[ic1] - fOffsetTAI) );
1619  delta += ( beamspill->GetT0Nano(i) - (8.*XCET1_coarse[ic1] + XCET1_frac[ic1] / 512.) );
1620 
1621  if( fXCETDebug ) std::cout << "XCET1 delta: " << delta << std::endl;
1622 
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;
1625  status_1.trigger = 1;
1626  status_1.timeStamp = delta;
1627  break;
1628  }
1629  }
1630  }
1631  beamspill->AddCKov0( status_1 );
1632 
1633  beam::CKov status_2;
1634  status_2.pressure = CKov2Pressure;
1635  status_2.timeStamp = -1.;
1636 
1637  if( !fetched_XCET2 ) status_2.trigger = -1;
1638  else{
1639 
1640  status_2.trigger = 0;
1641 
1642  for( size_t ic2 = 0; ic2 < XCET2_seconds.size(); ++ic2 ){
1643 
1644  double delta = 1.e9 * ( beamspill->GetT0Sec(i) - (XCET2_seconds[ic2] - fOffsetTAI) );
1645  delta += ( beamspill->GetT0Nano(i) - (8.*XCET2_coarse[ic2] + XCET2_frac[ic2] / 512.) );
1646 
1647  if( fXCETDebug ) std::cout << "XCET2 delta: " << delta << std::endl;
1648 
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;
1651  status_2.trigger = 1;
1652  status_2.timeStamp = delta;
1653  break;
1654  }
1655  }
1656  }
1657  beamspill->AddCKov1( status_2 );
1658 
1659  }
1660 
1661  if( fXCETDebug ){
1662  MF_LOG_INFO("BeamEvent") << "GeneralTriggers: " << beamspill->GetNT0() << std::endl;
1663  MF_LOG_INFO("BeamEvent") << "XCET1: " << beamspill->GetNCKov0() << std::endl;
1664  MF_LOG_INFO("BeamEvent") << "XCET2: " << beamspill->GetNCKov1() << std::endl;
1665 
1666  int nxcet1 = 0, nxcet2 = 0;
1667 
1668  for( size_t i = 0; i < beamspill->GetNT0(); ++i ){
1669  if( beamspill->GetCKov0Status(i) == 1 ) nxcet1++ ;
1670  if( beamspill->GetCKov1Status(i) == 1 ) nxcet2++ ;
1671  }
1672 
1673  if( nxcet1 != CKov1Counts ) MF_LOG_WARNING("BeamEvent") << "CKov1 counts differ. In spill: " << nxcet1 << " From counts: " << CKov1Counts << "\n";
1674  if( nxcet2 != CKov2Counts ) MF_LOG_WARNING("BeamEvent") << "CKov2 counts differ. In spill: " << nxcet2 << " From counts: " << CKov2Counts << "\n";
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;
1677  }
1678 
1679 }
1680 
1681 
1682 ////////////////////////
1683 //
1685 
1686  // Retrieve the number of counts in the BPF
1687  std::vector<double> counts;
1688 
1689  try{
1690  counts = FetchAndReport(time, fXBPFPrefix + name + ":countsRecords[]", bfp);
1691  }
1692  catch( std::exception const&){
1693  MF_LOG_WARNING("BeamEvent") << "Could not fetch " << fXBPFPrefix + name + ":countsRecords[]\n";
1694  return;
1695  }
1696 
1697  if( fPrintDebug )
1698  MF_LOG_INFO("BeamEvent") << "Counts: " << counts[1] << "\n";
1699 
1700  std::vector<double> data;
1701  try{
1702  data = FetchAndReport(time, fXBPFPrefix + name + ":eventsData[]", bfp);
1703  }
1704  catch( std::exception const&){
1705  MF_LOG_WARNING("BeamEvent") << "Could not fetch " << fXBPFPrefix + name + ":eventsData[]\n";
1706  return;
1707  }
1708 
1709  if( fPrintDebug ){
1710  // If the number of counts is larger than the number of general triggers
1711  // make note
1712  if(counts[1] != beamspill->GetNT0()){
1713  MF_LOG_WARNING("BeamEvent") << "WARNING MISMATCH " << counts[1] << " " << beamspill->GetNT0() << "\n";
1714  }
1715  }
1716 
1717  beam::FBM fbm(ID);
1718  /*
1719  beam::FBM fbm;
1720  fbm.ID = ID;
1721  fbm.fibers = {};
1722  fbm.glitch_mask = {};
1723  std::uninitialized_fill( std::begin(fbm.fiberData), std::end(fbm.fiberData), 0. );
1724  std::uninitialized_fill( std::begin(fbm.timeData), std::end(fbm.timeData), 0. );
1725  fbm.timeStamp = 0.;
1726  fbm.decoded = false;
1727  fbm.active = std::vector<short>();
1728  */
1729 
1730  //Use this just in case any are out of sync?
1731  //Shouldn't be, but just to be safe...
1732  //Helps cut down on time
1733  std::vector<size_t> leftOvers;
1734  for(size_t lo = 0; lo < beamspill->GetNT0(); ++lo){
1735  leftOvers.push_back(lo);
1736  }
1737 
1738  //std::cout.precision(20);
1739  for(size_t i = 0; i < counts[1]; ++i){
1740  for(int j = 0; j < 10; ++j){
1741 
1742  double theData = data[20*i + (2*j + 1)];
1743 
1744  //std::cout << "theData: " << theData << std::endl;
1745 
1746  if(j < 4)
1747  fbm.timeData[j] = theData;
1748  else
1749  fbm.fiberData[j - 4] = theData;
1750 
1751  }
1752 
1753  // Check the time data for corruption
1754  if(fbm.timeData[1] < .0000001)
1755  continue;
1756 
1757 
1758  for(std::vector<size_t>::iterator ip = leftOvers.begin(); ip != leftOvers.end(); ++ip){
1759 
1760  // Compute the time delta between the timeStamp and the T0, see if it's less than 500ns away
1761 
1762  double delta = (beamspill->GetT0Sec(*ip) - (fbm.timeData[3] - fOffsetTAI) );
1763  delta += 1.e-9*( beamspill->GetT0Nano(*ip) - 8.*fbm.timeData[2] );
1764 
1765  if( delta < -1.e-7 && delta > -10.e-7 ){
1766 
1767  beamspill->ReplaceFBMTrigger(name, fbm, *ip);
1768  leftOvers.erase(ip);
1769  break;
1770  }
1771  }
1772  }
1773 
1774  if( fPrintDebug ){
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){
1778  MF_LOG_WARNING("BeamEvent") << leftOvers[ip] << " ";
1779  }
1780  MF_LOG_WARNING("BeamEvent") << "\n";
1781  }
1782  }
1783 
1784  for(size_t i = 0; i < beamspill->GetNFBMTriggers(name); ++i){
1785  beamspill->DecodeFibers(name,i);
1786  }
1787 
1788  //Check the last 2 words for any duplicates
1789  beamspill->FixFiberGlitch(name);
1790  if( fPrintDebug ){
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;
1796  }
1797  }
1798  }
1799  }
1800 
1801 }
1802 
1803 ////////////////////////
1804 //
1806  for(size_t d = 0; d < fDevices.size(); ++d){
1808  parseGeneralXBPF(name, time, d);
1809  }
1810 }
1811 // END BeamEvent::parseXBFP
1812 ////////////////////////
1813 
1814 
1816 {
1818 
1819 
1820  if( fDebugMomentum ){
1821  fFullMomentum = tfs->make<TH1F>("FullMomentum", "", 150, 0, 15.);
1822  fCutMomentum = tfs->make<TH1F>("CutMomentum", "", 150, 0, 15.);
1823  }
1824 
1825  if( fSaveOutTree ){
1826  fOutTree = tfs->make<TTree>("tree", "lines");
1827  fOutTree->Branch("Time", &eventTime);
1828  fOutTree->Branch("RDTS", &RDTSTime);
1829  fOutTree->Branch("RDTSSec", &RDTSSec);
1830  fOutTree->Branch("RDTSNano", &RDTSNano);
1831  fOutTree->Branch("RDTSTrigger", &RDTSTrigger);
1832  fOutTree->Branch("Event", &eventNum);
1833  fOutTree->Branch("SpillStart", &SpillStart);
1834  fOutTree->Branch("SpillEnd", &SpillEnd);
1835  fOutTree->Branch("SpillOffset", &SpillOffset);
1836  fOutTree->Branch("ActiveTriggerTime", &ActiveTriggerTime);
1837  fOutTree->Branch("Run", &runNum);
1838  fOutTree->Branch("Subrun", &subRunNum);
1839  fOutTree->Branch("Pressure1", &CKov1Pressure);
1840  fOutTree->Branch("Counts1", &CKov1Counts);
1841  fOutTree->Branch("Pressure2", &CKov2Pressure);
1842  fOutTree->Branch("Counts2", &CKov2Counts);
1843  fOutTree->Branch("acqTime", &acqTime);
1844  fOutTree->Branch("acqStampMBPL", &acqStampMBPL);
1845  fOutTree->Branch("C1DB", &C1DB);
1846  fOutTree->Branch("C2DB", &C2DB);
1847  fOutTree->Branch("s11Sec", &s11Sec);
1848  fOutTree->Branch("s11Nano", &s11Nano);
1849  fOutTree->Branch("genTrigFracs", &genTrigFracs);
1850  fOutTree->Branch("genTrigCoarses", &genTrigCoarses);
1851  fOutTree->Branch("genTrigSecs", &genTrigSecs);
1852  fOutTree->Branch("matched", &fMatched);
1853 
1854  fOutTree->Branch("TOF", &fOutTOF );
1855  fOutTree->Branch("P", &fOutP );
1856  fOutTree->Branch("C0", &fOutC0 );
1857  fOutTree->Branch("C1", &fOutC1 );
1858  }
1859 
1860  if( fDebugTOFs ){
1861  fGenTrigTree = tfs->make<TTree>("GenTrig","");
1862  fGenTrigTree->Branch("coarse", &fGenTrigCoarse);
1863  fGenTrigTree->Branch("frac", &fGenTrigFrac);
1864  fGenTrigTree->Branch("sec", &fGenTrigSec);
1865 
1866  fXTOF1ATree = tfs->make<TTree>("TOF1A","");
1867  fXTOF1ATree->Branch("coarse", &fXTOF1ACoarse);
1868  fXTOF1ATree->Branch("frac", &fXTOF1AFrac);
1869  fXTOF1ATree->Branch("sec", &fXTOF1ASec);
1870 
1871  fXTOF1BTree = tfs->make<TTree>("TOF1B","");
1872  fXTOF1BTree->Branch("coarse", &fXTOF1BCoarse);
1873  fXTOF1BTree->Branch("frac", &fXTOF1BFrac);
1874  fXTOF1BTree->Branch("sec", &fXTOF1BSec);
1875 
1876  fXTOF2ATree = tfs->make<TTree>("TOF2A","");
1877  fXTOF2ATree->Branch("coarse", &fXTOF2ACoarse);
1878  fXTOF2ATree->Branch("frac", &fXTOF2AFrac);
1879  fXTOF2ATree->Branch("sec", &fXTOF2ASec);
1880  fXTOF2ATree->Branch("diff2A", &diff2A);
1881 
1882  fXTOF2BTree = tfs->make<TTree>("TOF2B","");
1883  fXTOF2BTree->Branch("coarse", &fXTOF2BCoarse);
1884  fXTOF2BTree->Branch("frac", &fXTOF2BFrac);
1885  fXTOF2BTree->Branch("sec", &fXTOF2BSec);
1886  fXTOF2BTree->Branch("diff2B", &diff2B);
1887  }
1888 
1889  //Tells IFBeam to print out debug statements
1890  ifbeam_ns::BeamFolder::_debug = fIFBeamDebug;
1891 
1893  if( fPrintDebug ){
1894  MF_LOG_INFO("BeamEvent") << "%%%%%%%%%% Got beam folder %%%%%%%%%%\n";
1895  MF_LOG_INFO("BeamEvent") << "%%%%%%%%%% Setting TimeWindow: " << fTimeWindow << " %%%%%%%%%%\n";
1896  }
1897 
1898  bfp->set_epsilon( fBFEpsilon );
1899 
1900  if( fPrintDebug )
1901  MF_LOG_INFO("BeamEvent") << "%%%%%%%%%% Set beam epislon " << fBFEpsilon << " %%%%%%%%%%\n";
1902 
1904 
1905  if( fPrintDebug ){
1906  MF_LOG_INFO("BeamEvent") << "%%%%%%%%%% Got beam folder %%%%%%%%%%\n";
1907  MF_LOG_INFO("BeamEvent") << "%%%%%%%%%% Setting TimeWindow: " << fTimeWindow << " %%%%%%%%%%\n";
1908  }
1909 
1910  bfp_xcet->set_epsilon( fXCETEpsilon );
1911 
1912  if( fPrintDebug )
1913  MF_LOG_INFO("BeamEvent") << "%%%%%%%%%% Set beam epislon " << fBFEpsilon << " %%%%%%%%%%\n";
1914 
1915 
1916  //Rotate the basis vectors of the FBMs
1918 }
1919 
1920 uint64_t proto::BeamEvent::joinHighLow(double high, double low){
1921 
1922  uint64_t low64 = (uint64_t)low;
1923 
1924  uint64_t high64 = (uint64_t)high;
1925 
1926  high64 = high64 << 32;
1927  uint64_t joined = high64 | low64;
1928 
1929  return joined;
1930 }
1931 
1932 TVector3 proto::BeamEvent::ConvertProfCoordinates(double x, double y, double z, double zOffset){
1933  double off = fNP04FrontZ - zOffset;
1934 
1935  TVector3 old(x,y,z);
1936 
1937  double newX = x*fBMBasisX.X() + y*fBMBasisY.X() /*+ (z-zOffset)*fBMBasisZ.X()*/ + off*fabs(fBMBasisZ.X());
1938  double newY = x*fBMBasisX.Y() + y*fBMBasisY.Y() /*+ (z-zOffset)*fBMBasisZ.Y()*/ + off*fabs(fBMBasisZ.Y());
1939  double newZ = x*fBMBasisX.Z() + y*fBMBasisY.Z() /*+ (z-zOffset) */ - off*fabs(fBMBasisZ.Z());
1940 
1941  newX += fBeamX*10.;
1942  newY += fBeamY*10.;
1943  newZ += fBeamZ*10.;
1944 
1945  TVector3 result(newX/10., newY/10., newZ/10.);
1946  return result;
1947 }
1948 
1953 
1954  rotated = true;
1955 }
1956 
1958  vec.RotateX( fRotateMonitorYZ * TMath::Pi()/180. );
1959  vec.RotateY( fRotateMonitorXZ * TMath::Pi()/180. );
1960 }
1961 
1962 void proto::BeamEvent::MakeTrack(size_t theTrigger){
1963 
1964  if( fPrintDebug )
1965  MF_LOG_INFO("BeamEvent") << "Making Track for time: " << beamspill->GetT0Sec(theTrigger) << " " << beamspill->GetT0Nano(theTrigger) << "\n";
1966 
1967  //Get the active fibers from the upstream tracking XBPF
1968  std::vector<short> firstUpstreamFibers = beamspill->GetActiveFibers(firstUpstreamName, theTrigger);
1969  std::vector<short> secondUpstreamFibers = beamspill->GetActiveFibers(secondUpstreamName, theTrigger);
1970 
1971  if( fPrintDebug ){
1972  MF_LOG_INFO("BeamEvent") << firstUpstreamName << " has " << firstUpstreamFibers.size() << " active fibers"<< "\n";
1973  for(size_t i = 0; i < firstUpstreamFibers.size(); ++i){
1974  MF_LOG_INFO("BeamEvent") << firstUpstreamFibers[i] << " ";
1975  }
1976  MF_LOG_INFO("BeamEvent") << "\n";
1977 
1978  MF_LOG_INFO("BeamEvent") << secondUpstreamName << " has " << secondUpstreamFibers.size() << " active fibers" << "\n";
1979  for(size_t i = 0; i < secondUpstreamFibers.size(); ++i){
1980  MF_LOG_INFO("BeamEvent") << secondUpstreamFibers[i] << " ";
1981  }
1982  MF_LOG_INFO("BeamEvent") << "\n";
1983  }
1984 
1985  std::array<short,192> firstUpstreamGlitches = beamspill->GetFBM( firstUpstreamName, theTrigger ).glitch_mask;
1986  std::array<short,192> secondUpstreamGlitches = beamspill->GetFBM( secondUpstreamName, theTrigger ).glitch_mask;
1987  MaskGlitches( firstUpstreamFibers, firstUpstreamGlitches );
1988  MaskGlitches( secondUpstreamFibers, secondUpstreamGlitches );
1989  if( fPrintDebug ) std::cout << firstUpstreamName << " After fix: " << firstUpstreamFibers.size() << std::endl;
1990  if( fPrintDebug ) std::cout << secondUpstreamName << " After fix: " << secondUpstreamFibers.size() << std::endl;
1991  //////////////////////////////////////////////
1992 
1993  //Get the active fibers from the downstream tracking XBPF
1994  std::vector<short> firstDownstreamFibers = beamspill->GetActiveFibers(firstDownstreamName, theTrigger);
1995  std::vector<short> secondDownstreamFibers = beamspill->GetActiveFibers(secondDownstreamName, theTrigger);
1996 
1997  if( fPrintDebug ){
1998  MF_LOG_INFO("BeamEvent") << firstDownstreamName << " has " << firstDownstreamFibers.size() << " active fibers" << "\n";
1999  for(size_t i = 0; i < firstDownstreamFibers.size(); ++i){
2000  MF_LOG_INFO("BeamEvent") << firstDownstreamFibers[i] << " ";
2001  }
2002  MF_LOG_INFO("BeamEvent") << "\n";
2003 
2004  MF_LOG_INFO("BeamEvent") << secondDownstreamName << " has " << secondDownstreamFibers.size() << " active fibers" << "\n";
2005  for(size_t i = 0; i < secondDownstreamFibers.size(); ++i){
2006  MF_LOG_INFO("BeamEvent") << secondDownstreamFibers[i] << " ";
2007  }
2008  MF_LOG_INFO("BeamEvent") << "\n";
2009  }
2010  std::array<short,192> firstDownstreamGlitches = beamspill->GetFBM( firstDownstreamName, theTrigger ).glitch_mask;
2011  std::array<short,192> secondDownstreamGlitches = beamspill->GetFBM( secondDownstreamName, theTrigger ).glitch_mask;
2012  MaskGlitches( firstDownstreamFibers, firstDownstreamGlitches );
2013  MaskGlitches( secondDownstreamFibers, secondDownstreamGlitches );
2014  if( fPrintDebug ) std::cout << firstDownstreamName << " After fix: " << firstDownstreamFibers.size() << std::endl;
2015  if( fPrintDebug ) std::cout << secondDownstreamName << " After fix: " << secondDownstreamFibers.size() << std::endl;
2016  //////////////////////////////////////////////
2017 
2018  if( (firstUpstreamFibers.size() < 1) || (secondUpstreamFibers.size() < 1) || (firstDownstreamFibers.size() < 1) || (secondDownstreamFibers.size() < 1) ){
2019  if( fPrintDebug )
2020  MF_LOG_INFO("BeamEvent") << "Warning, at least one empty Beam Profiler. Not making track" << "\n";
2021  return;
2022  }
2023 
2024  //We have the active Fibers, now go through them.
2025  //Skip the second of any adjacents
2026  std::vector< std::pair<size_t, size_t> > upstreamPairedFibers;
2027  std::vector< std::pair<size_t, size_t> > downstreamPairedFibers;
2028  std::string firstUpstreamType = fDeviceTypes[firstUpstreamName];
2029  std::string secondUpstreamType = fDeviceTypes[secondUpstreamName];
2030  std::string firstDownstreamType = fDeviceTypes[firstDownstreamName];
2031  std::string secondDownstreamType = fDeviceTypes[secondDownstreamName];
2032 
2033  std::vector< TVector3 > upstreamPositions;
2034  std::vector< TVector3 > downstreamPositions;
2035 
2036  //Pair the upstream fibers together
2037  if( fPrintDebug )
2038  MF_LOG_INFO("BeamEvent") << "Upstream" << "\n";
2039 
2040  for(size_t iF1 = 0; iF1 < firstUpstreamFibers.size(); ++iF1){
2041 
2042  size_t firstFiber = firstUpstreamFibers[iF1];
2043 
2044  for(size_t iF2 = 0; iF2 < secondUpstreamFibers.size(); ++iF2){
2045  size_t secondFiber = secondUpstreamFibers[iF2];
2046 
2047  if( fPrintDebug )
2048  MF_LOG_INFO("BeamEvent") << "Paired: " << firstFiber << " " << secondFiber << "\n";
2049 
2050  upstreamPairedFibers.push_back(std::make_pair(firstFiber, secondFiber));
2051 
2052  if (iF2 < secondUpstreamFibers.size() - 1){
2053  if (secondUpstreamFibers[iF2] == (secondUpstreamFibers[iF2 + 1] - 1)) ++iF2;
2054  }
2055  }
2056 
2057  if (iF1 < firstUpstreamFibers.size() - 1){
2058  if (firstUpstreamFibers[iF1] == (firstUpstreamFibers[iF1 + 1] - 1)) ++iF1;
2059  }
2060  }
2061 
2062  for(size_t iF = 0; iF < upstreamPairedFibers.size(); ++iF){
2063 
2064  std::pair<size_t,size_t> thePair = upstreamPairedFibers.at(iF);
2065 
2066  if(firstUpstreamType == "horiz" && secondUpstreamType == "vert"){
2067  double xPos = GetPosition(firstUpstreamName, thePair.first);
2068  double yPos = GetPosition(secondUpstreamName, thePair.second);
2069 
2070  TVector3 posInDet = ConvertProfCoordinates(xPos,yPos,0.,fFirstTrackingProfZ);
2071  upstreamPositions.push_back( posInDet );
2072 
2073  if( fPrintDebug ){
2074  MF_LOG_INFO("BeamEvent") << "normal " << xPos << " " << yPos << "\n";
2075  MF_LOG_INFO("BeamEvent") << posInDet.X() << " " << posInDet.Y() << " " << posInDet.Z() << "\n";
2076  }
2077 
2078  }
2079  else if(firstUpstreamType == "vert" && secondUpstreamType == "horiz"){
2080  double yPos = GetPosition(firstUpstreamName, thePair.first);
2081  double xPos = GetPosition(secondUpstreamName, thePair.second);
2082 
2083  TVector3 posInDet = ConvertProfCoordinates(xPos,yPos,0.,fFirstTrackingProfZ);
2084  upstreamPositions.push_back( posInDet );
2085 
2086  if( fPrintDebug ){
2087  MF_LOG_INFO("BeamEvent") << "normal " << xPos << " " << yPos << "\n";
2088  MF_LOG_INFO("BeamEvent") << posInDet.X() << " " << posInDet.Y() << " " << posInDet.Z() << "\n";
2089  }
2090  }
2091 
2092  }
2093 
2094  if( fPrintDebug )
2095  MF_LOG_INFO("BeamEvent") << "Downstream" << "\n";
2096 
2097  //Pair the downstream fibers together
2098  for(size_t iF1 = 0; iF1 < firstDownstreamFibers.size(); ++iF1){
2099 
2100  size_t firstFiber = firstDownstreamFibers[iF1];
2101 
2102  for(size_t iF2 = 0; iF2 < secondDownstreamFibers.size(); ++iF2){
2103  size_t secondFiber = secondDownstreamFibers[iF2];
2104 
2105  if( fPrintDebug )
2106  MF_LOG_INFO("BeamEvent") << "Paired: " << firstFiber << " " << secondFiber << "\n";
2107  downstreamPairedFibers.push_back(std::make_pair(firstFiber, secondFiber));
2108 
2109  if (iF2 < secondDownstreamFibers.size() - 1){
2110  if (secondDownstreamFibers[iF2] == (secondDownstreamFibers[iF2 + 1] - 1)) ++iF2;
2111  }
2112  }
2113 
2114  if (iF1 < firstDownstreamFibers.size() - 1){
2115  if (firstDownstreamFibers[iF1] == (firstDownstreamFibers[iF1 + 1] - 1)) ++iF1;
2116  }
2117  }
2118 
2119  for(size_t iF = 0; iF < downstreamPairedFibers.size(); ++iF){
2120 
2121  std::pair<size_t,size_t> thePair = downstreamPairedFibers.at(iF);
2122 
2123  if(firstDownstreamType == "horiz" && secondDownstreamType == "vert"){
2124  double xPos = GetPosition(firstDownstreamName, thePair.first);
2125  double yPos = GetPosition(secondDownstreamName, thePair.second);
2126 
2127  TVector3 posInDet = ConvertProfCoordinates(xPos,yPos,0.,fSecondTrackingProfZ);
2128  downstreamPositions.push_back( posInDet );
2129 
2130  if( fPrintDebug ){
2131  MF_LOG_INFO("BeamEvent") << "normal " << xPos << " " << yPos << "\n";
2132  MF_LOG_INFO("BeamEvent") << posInDet.X() << " " << posInDet.Y() << " " << posInDet.Z() << "\n";
2133  }
2134  }
2135  else if(firstDownstreamType == "vert" && secondDownstreamType == "horiz"){
2136  double yPos = GetPosition(firstDownstreamName, thePair.first);
2137  double xPos = GetPosition(secondDownstreamName, thePair.second);
2138 
2139  TVector3 posInDet = ConvertProfCoordinates(xPos,yPos,0.,fSecondTrackingProfZ);
2140  downstreamPositions.push_back( posInDet );
2141 
2142  if( fPrintDebug ){
2143  MF_LOG_INFO("BeamEvent") << "normal " << xPos << " " << yPos << "\n";
2144  MF_LOG_INFO("BeamEvent") << posInDet.X() << " " << posInDet.Y() << " " << posInDet.Z() << "\n";
2145  }
2146 
2147  }
2148 
2149  }
2150 
2151 
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));
2157 
2158  //Now project the last point to the TPC face
2159  thePoints.push_back( ProjectToTPC(thePoints[0],thePoints[1]) );
2160  std::vector<TVector3> theMomenta;
2161 
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() );
2165 
2168  recob::Track::Flags_t(thePoints.size()), false),
2170 
2171  beamevt->AddBeamTrack( tempTrack );
2172  }
2173  }
2174 
2175 }
2176 
2177 void proto::BeamEvent::MomentumSpec(size_t theTrigger){
2178 
2179  if( fPrintDebug )
2180  MF_LOG_INFO("BeamEvent") << "Doing momentum spectrometry for trigger " << beamspill->GetT0Sec(theTrigger) << " " << beamspill->GetT0Nano(theTrigger) << "\n";
2181 
2182  double LB = mag_P1*fabs(current[0]);
2183  double deltaI = fabs(current[0]) - mag_P4;
2184  if(deltaI>0) LB+= mag_P3*deltaI*deltaI;
2185 
2186  //Get the active fibers from the upstream tracking XBPF
2187  std::string firstBPROF1Type = fDeviceTypes[firstBPROF1];
2188  std::string secondBPROF1Type = fDeviceTypes[secondBPROF1];
2189  std::vector<short> BPROF1Fibers;
2190  std::string BPROF1Name;
2191 
2192  if (firstBPROF1Type == "horiz" && secondBPROF1Type == "vert"){
2193  BPROF1Fibers = beamspill->GetActiveFibers(firstBPROF1, theTrigger);
2194  BPROF1Name = firstBPROF1;
2195  }
2196  else if(secondBPROF1Type == "horiz" && firstBPROF1Type == "vert"){
2197  BPROF1Fibers = beamspill->GetActiveFibers(secondBPROF1, theTrigger);
2198  BPROF1Name = secondBPROF1;
2199  }
2200  else{
2201  MF_LOG_WARNING("BeamEvent") << "Error: type is not correct" << "\n";
2202  return;
2203  }
2204 
2205  if( fPrintDebug ){
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"
2214  << "\t" << beamspill->GetFBM( BPROF1Name, theTrigger ).fiberData[5] << std::endl;
2215 
2216  for(size_t i = 0; i < BPROF1Fibers.size(); ++i){
2217  MF_LOG_INFO("BeamEvent") << BPROF1Fibers[i] << " ";
2218  }
2219  MF_LOG_INFO("BeamEvent") << "\n";
2220  }
2221  //////////////////////////////////////////////
2222 
2223  //Remove the glitched fibers from BPROF1
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;
2226  MaskGlitches( BPROF1Fibers, BPROF1Glitches );
2227  if( fPrintDebug ) std::cout << "After fix: " << BPROF1Fibers.size() << std::endl;
2228 
2229 
2230 
2231  if( (BPROF1Fibers.size() < 1) ){
2232  if( fPrintDebug )
2233  MF_LOG_INFO("BeamEvent") << "Warning, at least one empty Beam Profiler. Not checking momentum" << "\n";
2234  return;
2235  }
2236  //We have the active Fibers, now go through them.
2237  //Skip the second of any adjacents
2238 
2239  //BPROF2////
2240  //
2241  std::vector<short> BPROF2Fibers = beamspill->GetActiveFibers(BPROF2, theTrigger);
2242 
2243  //Remove the glitched fibers from BPROF2
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;
2246  MaskGlitches( BPROF2Fibers, BPROF2Glitches );
2247  if( fPrintDebug ) std::cout << "After Fix " << BPROF2Fibers.size() << std::endl;
2248 
2249  if( (BPROF2Fibers.size() < 1) ){
2250  if( fPrintDebug )
2251  MF_LOG_INFO("BeamEvent") << "Warning, at least one empty Beam Profiler. Not checking momentum" << "\n";
2252  return;
2253  }
2254 
2255  if( fPrintDebug ){
2256  MF_LOG_INFO("BeamEvent") << BPROF2 << " has " << BPROF2Fibers.size() << " active fibers at time " << beamspill->GetFiberTime(BPROF2,theTrigger) << "\n";
2257  for(size_t i = 0; i < BPROF2Fibers.size(); ++i){
2258  MF_LOG_INFO("BeamEvent") << BPROF2Fibers[i] << " ";
2259  }
2260  MF_LOG_INFO("BeamEvent") << "\n";
2261  }
2262  ////////////
2263 
2264  //BPROF3////
2265  //
2266  std::vector<short> BPROF3Fibers = beamspill->GetActiveFibers(BPROF3, theTrigger);
2267 
2268  //Remove the glitched fibers from BPROF3
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;
2271  MaskGlitches( BPROF3Fibers, BPROF3Glitches );
2272  if( fPrintDebug ) std::cout << "After Fix " << BPROF3Fibers.size() << std::endl;
2273 
2274 
2275  if( (BPROF3Fibers.size() < 1) ){
2276  if( fPrintDebug )
2277  MF_LOG_INFO("BeamEvent") << "Warning, at least one empty Beam Profiler. Not checking momentum" << "\n";
2278  return;
2279  }
2280 
2281  if( fPrintDebug ){
2282  MF_LOG_INFO("BeamEvent") << BPROF3 << " has " << BPROF3Fibers.size() << " active fibers at time " << beamspill->GetFiberTime(BPROF3,theTrigger) << "\n";
2283  for(size_t i = 0; i < BPROF3Fibers.size(); ++i){
2284  MF_LOG_INFO("BeamEvent") << BPROF3Fibers[i] << " ";
2285  }
2286  MF_LOG_INFO("BeamEvent") << "\n";
2287  }
2288 
2289  if( fPrintDebug ){
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";
2294  }
2295 
2296  for(size_t i1 = 0; i1 < BPROF1Fibers.size(); ++i1){
2297  double x1,x2,x3;
2298 
2299  x1 = -1.*GetPosition(BPROF1Name, BPROF1Fibers[i1])/1.E3;
2300  if (i1 < BPROF1Fibers.size() - 1){
2301  if (BPROF1Fibers[i1] == (BPROF1Fibers[i1 + 1] - 1)){
2302  //Add .5 mm
2303  x1 += .0005;
2304  }
2305  }
2306 
2307  for(size_t i2 = 0; i2 < BPROF2Fibers.size(); ++i2){
2308  x2 = -1.*GetPosition(BPROF2, BPROF2Fibers[i2])/1.E3;
2309  if (i2 < BPROF2Fibers.size() - 1){
2310  if (BPROF2Fibers[i2] == (BPROF2Fibers[i2 + 1] - 1)){
2311  //Add .5 mm
2312  x2 += .0005;
2313  }
2314  }
2315 
2316  for(size_t i3 = 0; i3 < BPROF3Fibers.size(); ++i3){
2317 
2318  if( fPrintDebug )
2319  MF_LOG_INFO("BeamEvent") << "\t" << i1 << " " << i2 << " " << i3 << "\n";
2320 
2321  x3 = -1.*GetPosition(BPROF3, BPROF3Fibers[i3])/1.E3;
2322  if (i3 < BPROF3Fibers.size() - 1){
2323  if (BPROF3Fibers[i3] == (BPROF3Fibers[i3 + 1] - 1)){
2324  //Add .5 mm
2325  x3 += .0005;
2326  }
2327  }
2328 
2329 
2330  //Calibrate the positions
2331  //-1.*( FiberPos ) -> -1.*( FiberPos + ShiftDist )
2332  // = -1.*FiberPos - ShiftDist
2333  x1 = x1 - fBProf1Shift*1.e-3;
2334  x2 = x2 - fBProf2Shift*1.e-3;
2335  x3 = x3 - fBProf3Shift*1.e-3;
2336 
2337  double cosTheta_full = MomentumCosTheta(x1,x2,x3);
2338  double momentum_full = 299792458*LB/(1.E9 * acos(cosTheta_full));
2339  if( fDebugMomentum ) fFullMomentum->Fill(momentum_full);
2340 
2341  beamevt->AddRecoBeamMomentum(momentum_full);
2342 
2343  if (i3 < BPROF3Fibers.size() - 1){
2344  if (BPROF3Fibers[i3] == (BPROF3Fibers[i3 + 1] - 1)){
2345  //Skip the next
2346  ++i3;
2347  }
2348  }
2349 
2350  }
2351 
2352  if (i2 < BPROF2Fibers.size() - 1){
2353  if (BPROF2Fibers[i2] == (BPROF2Fibers[i2 + 1] - 1)){
2354  //Skip the next
2355  ++i2;
2356  }
2357  }
2358  }
2359 
2360  if (i1 < BPROF1Fibers.size() - 1){
2361  if (BPROF1Fibers[i1] == (BPROF1Fibers[i1 + 1] - 1)){
2362  //Skip the next
2363  ++i1;
2364  }
2365  }
2366  }
2367 }
2368 
2369 double proto::BeamEvent::MomentumCosTheta(double X1, double X2, double X3){
2370  double a = (X2*L3 - X3*L2)*cos(fBeamBend)/(L3-L2);
2371 
2372 
2373  double numTerm = (a - X1)*( (L3 - L2)*tan(fBeamBend) + (X3 - X2)*cos(fBeamBend) ) + L1*( L3 - L2 );
2374 
2375  double denomTerm1, denomTerm2, denom;
2376  denomTerm1 = sqrt( L1*L1 + (a - X1)*(a - X1) );
2377  denomTerm2 = sqrt( TMath::Power( ( (L3 - L2)*tan(fBeamBend) + (X3 - X2)*cos(fBeamBend) ),2)
2378  + TMath::Power( ( (L3 - L2) ),2) );
2379  denom = denomTerm1 * denomTerm2;
2380 
2381  double cosTheta = numTerm/denom;
2382  return cosTheta;
2383 }
2384 
2385 TVector3 proto::BeamEvent::ProjectToTPC(TVector3 firstPoint, TVector3 secondPoint){
2386  TVector3 dR = (secondPoint - firstPoint);
2387 
2388  double deltaZ = -1.*secondPoint.Z();
2389  double deltaX = deltaZ * (dR.X() / dR.Z());
2390  double deltaY = deltaZ * (dR.Y() / dR.Z());
2391 
2392  TVector3 lastPoint = secondPoint + TVector3(deltaX, deltaY, deltaZ);
2393  return lastPoint;
2394 }
2395 
2396 double proto::BeamEvent::GetPosition(std::string deviceName, int fiberIdx){
2397  if(fiberIdx > 192){ MF_LOG_WARNING("BeamEvent") << "Please select fiber in range [0,191]" << "\n"; return -1.;}
2398  double size = fFiberDimension[deviceName];
2399 
2400  //Define 0th fiber as farthest positive. Last fiber is farthest negative. Center is between 96 and 97
2401  double pos = size*(96 - fiberIdx) - size/2.;
2402  return pos;
2403 }
2404 
2405 void proto::BeamEvent::MaskGlitches( std::vector<short> & fibers, std::array<short,192> & glitches ){
2406  for( short i = 0; i < 192; ++i ){
2407  if( glitches[i] ){
2408  fibers.erase( std::find( fibers.begin(), fibers.end(), i ) );
2409  }
2410  }
2411 }
2412 
2413 
static QCString name
Definition: declinfo.cpp:673
std::string secondBPROF1
intermediate_table::iterator iterator
void GetSpillInfo(art::Event &)
std::string fXCETBundleName
double GetPosition(std::string, int)
void InitXBPFInfo(beam::ProtoDUNEBeamSpill *)
std::string fBundleName
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
EventNumber_t event() const
Definition: DataViewImpl.cc:85
const double & GetT0Nano() const
std::vector< std::string > fDevices
const int & GetTimingTrigger() const
void SetSpillStart(double theSpillStart)
static QCString result
std::vector< double > FetchAndReport(long long, std::string, std::unique_ptr< ifbeam_ns::BeamFolder > &)
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
void getS11Info(uint64_t)
unsigned int ID
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 > > SMatrixSym55
Definition: TrackingTypes.h:85
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:58
const double & GetCKov0Pressure() const
std::map< std::string, double > fFiberDimension
std::string fXCETPrefix
std::vector< double > diff2B
void SetTimingTrigger(int theTrigger)
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
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
Definition: Track.h:68
void MaskGlitches(std::vector< short > &, std::array< short, 192 > &)
#define MF_LOG_ERROR(category)
double GetPairedPosition(std::string, size_t)
STL namespace.
art::Handle< std::vector< raw::RDTimeStamp > > RDTimeStampHandle
uint8_t channel
Definition: CRTFragment.hh:201
beam::ProtoDUNEBeamEvent prev_beamevt
void MakeTrack(size_t)
void InitFBMs(std::vector< std::string >)
TVector3 ProjectToTPC(TVector3, TVector3)
beam::ProtoDUNEBeamEvent * beamevt
void parseXCETDB(uint64_t)
std::string fXBPFPrefix
BeamEvent & operator=(BeamEvent const &)=delete
long long int eventTime
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.
Definition: StdUtils.h:92
void SetCKov0(CKov theCKov)
double MomentumCosTheta(double, double, double)
T abs(T value)
const double e
void beginJob() override
Timestamp time() const
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
void SetMagnetCurrent(double theMagnetCurrent)
std::string fOutputLabel
std::vector< double > genTrigCoarses
A trajectory in space reconstructed from hits.
const double a
const double & GetTOF() const
const double & GetT0Sec() const
counts_as<> counts
Number of ADC counts, represented by signed short int.
Definition: electronics.h:116
def move(depos, offset)
Definition: depos.py:107
art::ServiceHandle< ifbeam_ns::IFBeam > ifb
void parseXTOF(uint64_t)
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::vector< double > diff2A
void SetRDTimestamp(long long theRDTimestamp)
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
p
Definition: test.py:223
uint64_t joinHighLow(double, double)
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
void SetCTBTimestamp(long long theCTBTimestamp)
General LArSoft Utilities.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:55
RunNumber_t run() const
Definition: DataViewImpl.cc:71
void AddBeamTrack(recob::Track theTrack)
#define MF_LOG_INFO(category)
std::string firstUpstreamName
void parseXBPF(uint64_t)
std::vector< double > current
beam::ProtoDUNEBeamSpill * beamspill
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
std::string fXTOFPrefix
void AddRecoBeamMomentum(double theMomentum)
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)
std::string firstBPROF1
void TimeIn(art::Event &, uint64_t)
const short & GetCKov0Status() const
list x
Definition: train.py:276
void SetT0(std::pair< double, double > theT0)
void parseGeneralXBPF(std::string, uint64_t, size_t)
void SetCKov1(CKov theCKov)
double fiberData[6]
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
#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
int bool
Definition: qglobal.h:345
void RotateMonitorVector(TVector3 &vec)
void SetTOFs(std::vector< double > theContent)
std::vector< double > genTrigFracs
double timeData[4]
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:
Definition: Track.h:49
std::vector< double > genTrigSecs
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
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)