Public Member Functions | Private Member Functions | Private Attributes | List of all members
protoana::ProtoDUNEBeamlineUtils Class Reference

#include <ProtoDUNEBeamlineUtils.h>

Public Member Functions

 ProtoDUNEBeamlineUtils (fhicl::ParameterSet const &pset)
 
 ~ProtoDUNEBeamlineUtils ()
 
void reconfigure (fhicl::ParameterSet const &pset)
 
const beam::ProtoDUNEBeamEvent GetBeamEvent (art::Event const &evt) const
 
void GetFibers (art::Event const &evt)
 
void GetCurrent (art::Event const &evt)
 
std::vector< recob::Track > MakeTracks (art::Event const &evt)
 
double GetPosition (short)
 
TVector3 ConvertMonitorCoordinates (double, double, double, double)
 
void BeamMonitorBasisVectors ()
 
void RotateMonitorVector (TVector3 &)
 
TVector3 ProjectToTPC (TVector3, TVector3)
 
std::vector< double > MomentumSpec (art::Event const &evt)
 
double MomentumCosTheta (double, double, double)
 
std::vector< int > GetPID (beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
 
std::vector< int > GetPID (art::Event const &evt, double nominal_momentum)
 
PossibleParticleCands GetPIDCandidates (beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
 
PossibleParticleCands GetPIDCandidates (art::Event const &beamevt, double nominal_momentum)
 
double ComputeMomentum (int pdg, double tof)
 
double ComputeTOF (int pdg, double momentum)
 
bool IsGoodBeamlineTrigger (art::Event const &evt) const
 
bool IsGoodBeamlineTrigger (beam::ProtoDUNEBeamEvent const &beamEvent) const
 
bool HasPerfectBeamMomentum (art::Event const &evt) const
 
bool HasPerfectBeamMomentum (beam::ProtoDUNEBeamEvent const &beamEvent) const
 
std::vector< double > GetBeamlineMass (art::Event const &evt) const
 
std::vector< double > GetBeamlineMassSquared (art::Event const &evt) const
 
const std::tuple< double, double, int, int > GetBeamlineVars (art::Event const &evt) const
 
const std::tuple< double, double, int, int, int, double, double, int, int, bool > GetBeamlineVarsAndStatus (art::Event const &evt) const
 

Private Member Functions

std::vector< int > GetPID_CERNCalib (beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
 
PossibleParticleCands GetPIDCandidates_CERNCalib (beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
 

Private Attributes

art::InputTag fBeamEventTag
 
bool fUseCERNCalibSelection
 
std::map< std::string, std::vector< short > > ActiveFibers
 
double Current
 
std::vector< std::string > AllDevices
 
std::string HorizUpstream = "XBPF022707"
 
std::string VertUpstream = "XBPF022708"
 
std::string HorizDownstream = "XBPF022716"
 
std::string VertDownstream = "XBPF022717"
 
std::string BProf1 = "XBPF022697"
 
std::string BProf2 = "XBPF022701"
 
std::string BProf3 = "XBPF022702"
 
TVector3 MonitorBasisX
 
TVector3 MonitorBasisY
 
TVector3 MonitorBasisZ
 
bool rotated = false
 
double fRotateMonitorXZ
 
double fRotateMonitorYZ
 
double fRotateMonitorYX
 
double fFirstTrackingProfZ
 
double fSecondTrackingProfZ
 
double fNP04FrontZ
 
double fBeamX
 
double fBeamY
 
double fBeamZ
 
double fBeamBend
 
double L1
 
double L2
 
double L3
 
double fBProf1Shift
 
double fBProf2Shift
 
double fBProf3Shift
 
double mag_P1 = 5.82044830e-3
 
double mag_P3 = -4.68880000e-6
 
double mag_P4 = 324.573967
 
std::map< int, double > particle_mass
 
double c = 299792458.
 
const double fTOFDist = 28.575
 
float fMomentumScaleFactor
 
float fMomentumOffset
 
float fTOFScaleFactor
 
float fTOFOffset
 

Detailed Description

Definition at line 76 of file ProtoDUNEBeamlineUtils.h.

Constructor & Destructor Documentation

protoana::ProtoDUNEBeamlineUtils::ProtoDUNEBeamlineUtils ( fhicl::ParameterSet const &  pset)

Definition at line 10 of file ProtoDUNEBeamlineUtils.cxx.

10  {
11  this->reconfigure(p);
12 }
test.p
p
Definition: test.py:223
protoana::ProtoDUNEBeamlineUtils::reconfigure
void reconfigure(fhicl::ParameterSet const &pset)
Definition: ProtoDUNEBeamlineUtils.cxx:205
protoana::ProtoDUNEBeamlineUtils::~ProtoDUNEBeamlineUtils ( )

Definition at line 14 of file ProtoDUNEBeamlineUtils.cxx.

14  {
15 
16 }

Member Function Documentation

void protoana::ProtoDUNEBeamlineUtils::BeamMonitorBasisVectors ( )

Definition at line 267 of file ProtoDUNEBeamlineUtils.cxx.

267  {
268  MonitorBasisX = TVector3( 1., 0., 0. );
269  MonitorBasisY = TVector3( 0., 1., 0. );
270  MonitorBasisZ = TVector3( 0., 0., 1. );
271 
272  RotateMonitorVector(MonitorBasisX);
273  RotateMonitorVector(MonitorBasisY);
274  RotateMonitorVector(MonitorBasisZ);
275 
276 
277  rotated = true;
278 }
double protoana::ProtoDUNEBeamlineUtils::ComputeMomentum ( int  pdg,
double  tof 
)

Compute what the beamline momentum SHOULD BE given a certain pdg and beamline TOF

Definition at line 632 of file ProtoDUNEBeamlineUtils.cxx.

632  {
633  if( tof < .0000001 ){
634  std::cout << "Low TOF" << std::endl;
635  return -1.;
636  }
637 
638  if( particle_mass.find( pdg ) == particle_mass.end() ){
639  std::cout << "PDG " << pdg << " not found" << std::endl;
640  return -1.;
641  }
642 
643  double m = particle_mass[pdg];
644  return m * sqrt(1. / ( 1. - fTOFDist*fTOFDist / ( 1.e-18*tof*tof*c*c ) ) - 1. );
645 
646 }
e
const double e
Definition: gUpMuFluxGen.cxx:165
cache_state.m
m
Definition: cache_state.py:415
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
double protoana::ProtoDUNEBeamlineUtils::ComputeTOF ( int  pdg,
double  momentum 
)

Compute what the TOF SHOULD BE given a certain pdg and beamline momentum

Definition at line 614 of file ProtoDUNEBeamlineUtils.cxx.

614  {
615  if( momentum < .0000001 ){
616  std::cout << "Low Momentum" << std::endl;
617  return -1.;
618  }
619 
620  if( particle_mass.find( pdg ) == particle_mass.end() ){
621  std::cout << "PDG " << pdg << " not found" << std::endl;
622  return -1.;
623  }
624 
625  double m = particle_mass[pdg];
626  double tof = fTOFDist / c;
627  tof = tof / sqrt( 1. - m*m / (momentum*momentum + m*m) );
628 
629  return tof * 1.e9;
630 }
cache_state.m
m
Definition: cache_state.py:415
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
reco_momentum_tuples.momentum
def momentum(x1, x2, x3, scale=1.)
Definition: reco_momentum_tuples.py:198
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
TVector3 protoana::ProtoDUNEBeamlineUtils::ConvertMonitorCoordinates ( double  x,
double  y,
double  z,
double  zOffset 
)

Definition at line 247 of file ProtoDUNEBeamlineUtils.cxx.

247  {
248 
249  if( !rotated ) BeamMonitorBasisVectors();
250 
251  double off = fNP04FrontZ - zOffset;
252 
253  TVector3 old(x,y,z);
254 
255  double newX = x*MonitorBasisX.X() + y*MonitorBasisY.X() + off*fabs(MonitorBasisZ.X());
256  double newY = x*MonitorBasisX.Y() + y*MonitorBasisY.Y() + off*fabs(MonitorBasisZ.Y());
257  double newZ = x*MonitorBasisX.Z() + y*MonitorBasisY.Z() - off*fabs(MonitorBasisZ.Z());
258 
259  newX += fBeamX*10.;
260  newY += fBeamY*10.;
261  newZ += fBeamZ*10.;
262 
263  TVector3 result(newX/10., newY/10., newZ/10.);
264  return result;
265 }
result
static QCString result
Definition: fortranscanner.cpp:56614
train.x
list x
Definition: train.py:276
const beam::ProtoDUNEBeamEvent protoana::ProtoDUNEBeamlineUtils::GetBeamEvent ( art::Event const &  evt) const

Definition at line 18 of file ProtoDUNEBeamlineUtils.cxx.

18  {
19  std::vector< art::Ptr< beam::ProtoDUNEBeamEvent > > beamVec;
20  auto beamHandle = evt.getValidHandle< std::vector< beam::ProtoDUNEBeamEvent> >(fBeamEventTag);
21  if( beamHandle.isValid() ){
22  art::fill_ptr_vector( beamVec, beamHandle );
23  }
24  const beam::ProtoDUNEBeamEvent & beamEvent = *(beamVec.at(0));
25  return beamEvent;
26 }
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
std::vector< double > protoana::ProtoDUNEBeamlineUtils::GetBeamlineMass ( art::Event const &  evt) const

Uses the beamline momentum and time of flight measurements to estimate the beamline particle mass in GeV/c^2

Returns a list of all possible cominations (there can be multiple momenta).

Definition at line 695 of file ProtoDUNEBeamlineUtils.cxx.

695  {
696  std::vector<double> result;
697 
698  const auto & massSquareds = GetBeamlineMassSquared(evt);
699  for (const auto & massSquared: massSquareds)
700  {
701  const auto & mass = std::sqrt(massSquared);
702  result.push_back(mass);
703  }
704  return result;
705 }
result
static QCString result
Definition: fortranscanner.cpp:56614
protoana::ProtoDUNEBeamlineUtils::GetBeamlineMassSquared
std::vector< double > GetBeamlineMassSquared(art::Event const &evt) const
Definition: ProtoDUNEBeamlineUtils.cxx:707
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
std::vector< double > protoana::ProtoDUNEBeamlineUtils::GetBeamlineMassSquared ( art::Event const &  evt) const

Uses the beamline momentum and time of flight measurements to estimate the beamline particle mass^2 in (GeV/c^2)^2

This can be more useful than the mass since measurement fluctuations can make it negative.

Returns a list of all possible cominations (there can be multiple momenta).

Definition at line 707 of file ProtoDUNEBeamlineUtils.cxx.

707  {
708  std::vector<double> tofs;
709  std::vector<double> momenta;
710  std::vector<double> result;
711 
712  std::vector<art::Ptr<beam::ProtoDUNEBeamEvent>> beamVec;
713  auto beamHand = evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>(fBeamEventTag);
714  if(beamHand.isValid())
715  {
716  art::fill_ptr_vector(beamVec, beamHand);
717  }
718 
719  for(size_t iBeamEvent=0; iBeamEvent < beamVec.size(); iBeamEvent++)
720  {
721  const beam::ProtoDUNEBeamEvent& beamEvent = *(beamVec.at(iBeamEvent));
722  tofs.push_back(beamEvent.GetTOF());
723 
724  const std::vector<double> & beamMomenta = beamEvent.GetRecoBeamMomenta();
725  for(size_t iMom=0; iMom < beamMomenta.size(); iMom++)
726  {
727  momenta.push_back(beamEvent.GetRecoBeamMomentum(iMom));
728  }
729  }
730  for(const auto & tof : tofs)
731  {
732  for(const auto & momentum : momenta)
733  {
734  const double massSquared = std::pow(momentum*fMomentumScaleFactor-fMomentumOffset,2)
735  * ((tof*fTOFScaleFactor-fTOFOffset)/(fTOFDist/c*1e9) - 1.);
736  result.push_back(massSquared);
737  }
738  }
739  return result;
740 }
result
static QCString result
Definition: fortranscanner.cpp:56614
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
beam::ProtoDUNEBeamEvent::GetRecoBeamMomenta
const std::vector< double > & GetRecoBeamMomenta() const
Definition: ProtoDUNEBeamEvent.h:91
beam::ProtoDUNEBeamEvent::GetTOF
const double & GetTOF() const
Definition: ProtoDUNEBeamEvent.h:66
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
reco_momentum_tuples.momentum
def momentum(x1, x2, x3, scale=1.)
Definition: reco_momentum_tuples.py:198
beam::ProtoDUNEBeamEvent::GetRecoBeamMomentum
const double & GetRecoBeamMomentum(size_t i) const
Definition: ProtoDUNEBeamEvent.h:93
const std::tuple< double, double, int, int > protoana::ProtoDUNEBeamlineUtils::GetBeamlineVars ( art::Event const &  evt) const

Get reconstructed beamline momentum (in GeV/c), tof (in ns), and flags for if the ckov's fired. Will be < 0 if invalid.

C++ structured binding, so call like: const auto [momentum, tof, ckov0, ckov1] = dataUtils.GetBeamlineInformation(e);

then you have the normal float momentum, int ckov0 variables, etc. in the current scope.

Definition at line 742 of file ProtoDUNEBeamlineUtils.cxx.

742  {
743 
744  double momentum = -99999.;
745  double tof = -99999.;
746  int ckov0 = -99999.;
747  int ckov1 = -99999.;
748 
749  std::vector<art::Ptr<beam::ProtoDUNEBeamEvent>> beamVec;
750  auto beamHand = evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>(fBeamEventTag);
751  if(beamHand.isValid())
752  {
753  art::fill_ptr_vector(beamVec, beamHand);
754  }
755  for(size_t iBeamEvent=0; iBeamEvent < beamVec.size(); iBeamEvent++)
756  {
757  const beam::ProtoDUNEBeamEvent& beamEvent = *(beamVec.at(iBeamEvent));
758  tof = beamEvent.GetTOF();
759  ckov0 = beamEvent.GetCKov0Status();
760  ckov1 = beamEvent.GetCKov1Status();
761 
762  const std::vector<double> & beamMomenta = beamEvent.GetRecoBeamMomenta();
763  if (beamMomenta.size() > 0)
764  {
765  momentum = beamEvent.GetRecoBeamMomentum(0);
766  }
767  }
768  return std::make_tuple(momentum, tof, ckov0, ckov1);
769 }
beam::ProtoDUNEBeamEvent::GetCKov1Status
const short & GetCKov1Status() const
Definition: ProtoDUNEBeamEvent.h:38
beam::ProtoDUNEBeamEvent::GetRecoBeamMomenta
const std::vector< double > & GetRecoBeamMomenta() const
Definition: ProtoDUNEBeamEvent.h:91
beam::ProtoDUNEBeamEvent::GetTOF
const double & GetTOF() const
Definition: ProtoDUNEBeamEvent.h:66
beam::ProtoDUNEBeamEvent::GetCKov0Status
const short & GetCKov0Status() const
Definition: ProtoDUNEBeamEvent.h:37
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
reco_momentum_tuples.momentum
def momentum(x1, x2, x3, scale=1.)
Definition: reco_momentum_tuples.py:198
beam::ProtoDUNEBeamEvent::GetRecoBeamMomentum
const double & GetRecoBeamMomentum(size_t i) const
Definition: ProtoDUNEBeamEvent.h:93
const std::tuple< double, double, int, int, int, double, double, int, int, bool > protoana::ProtoDUNEBeamlineUtils::GetBeamlineVarsAndStatus ( art::Event const &  evt) const

Get reconstructed beamline momentum (in GeV/c), tof (in ns), flags for the tofChannel and ckov's. Also the timing trigger (12 means beam), the beam instrumentation trigger (1 means beam trigger), and whether the timing and beam instrumentation triggers are matched.

All values will be < 0 if invalid.

C++ structured binding, so call like: const auto [momentum, tof, tofChannel,ckov0,ckov1,ckov0Pressure,ckov1Pressure,timingTrigger,BITrigger,areBIAndTimingMatched] = dataUtils.GetBeamlineInformation(e);

then you have the normal float momentum, int ckov0 variables, etc. in the current scope. You can use (void) variable; lines to get rid of unused var warnings.

Definition at line 771 of file ProtoDUNEBeamlineUtils.cxx.

771  {
772 
773  double momentum = -99999.;
774  double tof = -99999.;
775  int tofChannel = -99999.;
776  int ckov0 = -99999.;
777  int ckov1 = -99999.;
778  double ckov0Pressure = -99999.;
779  double ckov1Pressure = -99999.;
780  int timingTrigger = -99999.;
781  int BITrigger = -99999.;
782  bool areBIAndTimingMatched = false;
783 
784  std::vector<art::Ptr<beam::ProtoDUNEBeamEvent>> beamVec;
785  auto beamHand = evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>(fBeamEventTag);
786  if(beamHand.isValid())
787  {
788  art::fill_ptr_vector(beamVec, beamHand);
789  }
790  for(size_t iBeamEvent=0; iBeamEvent < beamVec.size(); iBeamEvent++)
791  {
792  const beam::ProtoDUNEBeamEvent& beamEvent = *(beamVec.at(iBeamEvent));
793  tof = beamEvent.GetTOF();
794  tofChannel = beamEvent.GetTOFChan();
795  ckov0 = beamEvent.GetCKov0Status();
796  ckov1 = beamEvent.GetCKov1Status();
797  ckov0Pressure = beamEvent.GetCKov0Pressure();
798  ckov1Pressure = beamEvent.GetCKov1Pressure();
799  timingTrigger = beamEvent.GetTimingTrigger();
800  BITrigger = beamEvent.GetBITrigger();
801  areBIAndTimingMatched = beamEvent.CheckIsMatched();
802 
803  const std::vector<double> & beamMomenta = beamEvent.GetRecoBeamMomenta();
804  if (beamMomenta.size() > 0)
805  {
806  momentum = beamEvent.GetRecoBeamMomentum(0);
807  }
808  break;
809  }
810  return std::make_tuple(momentum, tof, tofChannel, ckov0, ckov1, ckov0Pressure, ckov1Pressure, timingTrigger, BITrigger, areBIAndTimingMatched);
811 }
beam::ProtoDUNEBeamEvent::GetTimingTrigger
const int & GetTimingTrigger() const
Definition: ProtoDUNEBeamEvent.h:100
beam::ProtoDUNEBeamEvent::GetCKov0Pressure
const double & GetCKov0Pressure() const
Definition: ProtoDUNEBeamEvent.h:41
beam::ProtoDUNEBeamEvent::GetCKov1Status
const short & GetCKov1Status() const
Definition: ProtoDUNEBeamEvent.h:38
beam::ProtoDUNEBeamEvent::GetRecoBeamMomenta
const std::vector< double > & GetRecoBeamMomenta() const
Definition: ProtoDUNEBeamEvent.h:91
beam::ProtoDUNEBeamEvent::GetTOF
const double & GetTOF() const
Definition: ProtoDUNEBeamEvent.h:66
beam::ProtoDUNEBeamEvent::GetCKov0Status
const short & GetCKov0Status() const
Definition: ProtoDUNEBeamEvent.h:37
beam::ProtoDUNEBeamEvent::GetBITrigger
const int & GetBITrigger() const
Definition: ProtoDUNEBeamEvent.h:97
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
reco_momentum_tuples.momentum
def momentum(x1, x2, x3, scale=1.)
Definition: reco_momentum_tuples.py:198
beam::ProtoDUNEBeamEvent::GetRecoBeamMomentum
const double & GetRecoBeamMomentum(size_t i) const
Definition: ProtoDUNEBeamEvent.h:93
beam::ProtoDUNEBeamEvent::CheckIsMatched
const bool & CheckIsMatched() const
Definition: ProtoDUNEBeamEvent.h:86
beam::ProtoDUNEBeamEvent::GetTOFChan
const int & GetTOFChan() const
Definition: ProtoDUNEBeamEvent.h:67
beam::ProtoDUNEBeamEvent::GetCKov1Pressure
const double & GetCKov1Pressure() const
Definition: ProtoDUNEBeamEvent.h:42
void protoana::ProtoDUNEBeamlineUtils::GetCurrent ( art::Event const &  evt)

Definition at line 27 of file ProtoDUNEBeamlineUtils.cxx.

27  {
28 
29  auto beamHandle = evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>(fBeamEventTag);
30 
31  std::vector<art::Ptr<beam::ProtoDUNEBeamEvent>> beamVec;
32  if( beamHandle.isValid()){
33  art::fill_ptr_vector(beamVec, beamHandle);
34  }
35 
36  //Should just have one
37  const beam::ProtoDUNEBeamEvent & beamEvent = *(beamVec.at(0));
38 
39  Current = beamEvent.GetMagnetCurrent();
40 }
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
beam::ProtoDUNEBeamEvent::GetMagnetCurrent
const double & GetMagnetCurrent() const
Definition: ProtoDUNEBeamEvent.h:115
void protoana::ProtoDUNEBeamlineUtils::GetFibers ( art::Event const &  evt)

Definition at line 42 of file ProtoDUNEBeamlineUtils.cxx.

42  {
43 
44  auto beamHandle = evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>(fBeamEventTag);
45 
46  std::vector<art::Ptr<beam::ProtoDUNEBeamEvent>> beamVec;
47  if( beamHandle.isValid()){
48  art::fill_ptr_vector(beamVec, beamHandle);
49  }
50 
51  //Should just have one
52  const beam::ProtoDUNEBeamEvent & beamEvent = *(beamVec.at(0));
53 
54  for( size_t i = 0; i < AllDevices.size(); ++i ){
55  std::string monitor = AllDevices[i];
56  ActiveFibers[ monitor ] = beamEvent.GetActiveFibers( monitor );
57 
58 
59  if( ActiveFibers.at(monitor).size() > 0){
60  std::cout << monitor << " has active fibers: " << std::endl;
61  for( size_t j = 0; j < ActiveFibers.at(monitor).size(); ++j ){
62  std::cout << ActiveFibers.at(monitor).at(j) << " ";
63  }
64  std::cout << std::endl << std::endl;
65  }
66  else{
67  std::cout << monitor << " has no active fibers" << std::endl;
68  }
69  }
70 }
protoana::ProtoDUNEBeamlineUtils::AllDevices
std::vector< std::string > AllDevices
Definition: ProtoDUNEBeamlineUtils.h:235
string
std::string string
Definition: nybbler.cc:12
beam::ProtoDUNEBeamEvent::GetActiveFibers
const std::vector< short > & GetActiveFibers(std::string) const
Definition: ProtoDUNEBeamEvent.h:257
protoana::ProtoDUNEBeamlineUtils::ActiveFibers
std::map< std::string, std::vector< short > > ActiveFibers
Definition: ProtoDUNEBeamlineUtils.h:227
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
std::vector< int > protoana::ProtoDUNEBeamlineUtils::GetPID ( beam::ProtoDUNEBeamEvent const &  beamevt,
double  nominal_momentum 
)

Get the particle ID from beamline instrumentation info (returns PDG IDs)

Make sure to set fcl param UseCERNCalibSelection=true if the beam reco is done after ~v08_07_00 otherwise if using an old reco'd BeamEvent (e.g. you are using the first central reco) set it to false

Definition at line 459 of file ProtoDUNEBeamlineUtils.cxx.

459  {
460  const auto& thePIDCands = GetPIDCandidates(beamevt, nominal_momentum);
461  std::vector< int > thePIDs = thePIDCands.getPDGCodes();
462  return thePIDs;
463 }
protoana::ProtoDUNEBeamlineUtils::GetPIDCandidates
PossibleParticleCands GetPIDCandidates(beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
Definition: ProtoDUNEBeamlineUtils.cxx:445
std::vector< int > protoana::ProtoDUNEBeamlineUtils::GetPID ( art::Event const &  evt,
double  nominal_momentum 
)

Definition at line 465 of file ProtoDUNEBeamlineUtils.cxx.

465  {
466  auto beamHand = evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>(fBeamEventTag);
467  for(size_t iBeamEvent=0; iBeamEvent < beamHand->size(); iBeamEvent++)
468  {
469  return GetPID(beamHand->at(iBeamEvent),nominal_momentum);
470  }
471  return std::vector<int>();
472 }
protoana::ProtoDUNEBeamlineUtils::GetPID
std::vector< int > GetPID(beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
Definition: ProtoDUNEBeamlineUtils.cxx:459
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
std::vector< int > protoana::ProtoDUNEBeamlineUtils::GetPID_CERNCalib ( beam::ProtoDUNEBeamEvent const &  beamevt,
double  nominal_momentum 
)
private

Get the particle ID from beamline instrumentation info (this uses the official CERN cuts and assumes CERN calibrations)

Definition at line 608 of file ProtoDUNEBeamlineUtils.cxx.

608  {
609  const auto& thePIDCands = GetPIDCandidates_CERNCalib(beamevt, nominal_momentum);
610  std::vector< int > thePIDs = thePIDCands.getPDGCodes();
611  return thePIDs;
612 }
protoana::ProtoDUNEBeamlineUtils::GetPIDCandidates_CERNCalib
PossibleParticleCands GetPIDCandidates_CERNCalib(beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
Definition: ProtoDUNEBeamlineUtils.cxx:474
protoana::PossibleParticleCands protoana::ProtoDUNEBeamlineUtils::GetPIDCandidates ( beam::ProtoDUNEBeamEvent const &  beamevt,
double  nominal_momentum 
)

Get the particle ID from beamline instrumentation info

Make sure to set fcl param UseCERNCalibSelection=true if the beam reco is done after ~v08_07_00 otherwise if using an old reco'd BeamEvent (e.g. you are using the first central reco) set it to false

Definition at line 445 of file ProtoDUNEBeamlineUtils.cxx.

445  {
446  return GetPIDCandidates_CERNCalib(beamevt,nominal_momentum);
447 }
protoana::ProtoDUNEBeamlineUtils::GetPIDCandidates_CERNCalib
PossibleParticleCands GetPIDCandidates_CERNCalib(beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
Definition: ProtoDUNEBeamlineUtils.cxx:474
protoana::PossibleParticleCands protoana::ProtoDUNEBeamlineUtils::GetPIDCandidates ( art::Event const &  beamevt,
double  nominal_momentum 
)

Definition at line 449 of file ProtoDUNEBeamlineUtils.cxx.

449  {
450  auto beamHand = evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>(fBeamEventTag);
451  for(size_t iBeamEvent=0; iBeamEvent < beamHand->size(); iBeamEvent++)
452  {
453  return GetPIDCandidates(beamHand->at(iBeamEvent),nominal_momentum);
454  }
455  return protoana::PossibleParticleCands();
456 }
protoana::ProtoDUNEBeamlineUtils::GetPIDCandidates
PossibleParticleCands GetPIDCandidates(beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
Definition: ProtoDUNEBeamlineUtils.cxx:445
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
protoana::PossibleParticleCands protoana::ProtoDUNEBeamlineUtils::GetPIDCandidates_CERNCalib ( beam::ProtoDUNEBeamEvent const &  beamevt,
double  nominal_momentum 
)
private

Get the particle ID from beamline instrumentation info (this uses the official CERN cuts and assumes CERN calibrations) if fUseCERNCalibSelection is set, otherwise uses old values.

Definition at line 474 of file ProtoDUNEBeamlineUtils.cxx.

474  {
475 
476  PossibleParticleCands candidates;
477 
478  //Check if momentum is in valid set
479  std::vector< double > valid_momenta = {1., 2., 3., 6., 7.};
480  if( std::find(valid_momenta.begin(), valid_momenta.end(), nominal_momentum) == valid_momenta.end() ){
481  std::cout << "Reference momentum " << nominal_momentum << " not valid" << std::endl;
482  return candidates;
483  }
484 
485  //Get the high/low pressure Cerenkov info
486  //int high_pressure_status, low_pressure_status;
487  int high_pressure_status = beamevt.GetCKov0Status();
488  int low_pressure_status = beamevt.GetCKov1Status();
489 
490  //std::cout << "Pressures: " << beamevt.GetCKov0Pressure() << " " << beamevt.GetCKov1Pressure() << std::endl;
491  //if( beamevt.GetCKov0Pressure() < beamevt.GetCKov1Pressure() ){
492  // high_pressure_status = beamevt.GetCKov1Status();
493  // low_pressure_status = beamevt.GetCKov0Status();
494  //}
495  //else{
496  // high_pressure_status = beamevt.GetCKov0Status();
497  // low_pressure_status = beamevt.GetCKov1Status();
498  //}
499 
500 
501 
502 
503  if( nominal_momentum == 1. ){
504  if( beamevt.GetTOFChan() == -1 ){
505  std::cout << "TOF invalid" << std::endl;
506  return candidates;
507  }
508  if( low_pressure_status == -1 ){
509  std::cout << "High pressure status invalid" << std::endl;
510  return candidates;
511  }
512 
513  const double & tof = beamevt.GetTOF();
514  if (
515  ((fUseCERNCalibSelection && tof < 105.)
516  || (!fUseCERNCalibSelection && tof < 170.))
517  && low_pressure_status == 1
518  ) {
519  candidates.electron = true;
520  }
521  else if (
522  ((fUseCERNCalibSelection && tof < 110.)
523  || (!fUseCERNCalibSelection && tof < 170.))
524  && low_pressure_status == 0 ){
525  candidates.muon = true;
526  candidates.pion = true;
527  }
528  else if (
529  ((fUseCERNCalibSelection && tof > 110. && tof < 160.)
530  || (!fUseCERNCalibSelection && tof > 170.))
531  && low_pressure_status == 0 ) {
532  candidates.proton = true;
533  }
534  }
535  else if( nominal_momentum == 2. ){
536  if( beamevt.GetTOFChan() == -1 ){
537  std::cout << "TOF invalid" << std::endl;
538  return candidates;
539  }
540  if( low_pressure_status == -1 ){
541  std::cout << "High pressure Cerenkov status invalid" << std::endl;
542  return candidates;
543  }
544 
545  const double & tof = beamevt.GetTOF();
546  if (
547  ((fUseCERNCalibSelection && tof < 105.)
548  || (!fUseCERNCalibSelection && tof < 160.))
549  && low_pressure_status == 1
550  ) {
551  candidates.electron = true;
552  }
553  else if (
554  ((fUseCERNCalibSelection && tof < 103.)
555  || (!fUseCERNCalibSelection && tof < 160.))
556  && low_pressure_status == 0 ){
557  candidates.muon = true;
558  candidates.pion = true;
559  }
560  else if (
561  ((fUseCERNCalibSelection && tof > 103. && tof < 160.)
562  || (!fUseCERNCalibSelection && tof > 160.))
563  && low_pressure_status == 0 ) {
564  candidates.proton = true;
565  }
566  }
567  else if( nominal_momentum == 3. ){
568  if( high_pressure_status == -1 || low_pressure_status == -1 ){
569  std::cout << "At least one Cerenkov status invalid " << std::endl;
570  std::cout << "High: " << high_pressure_status << " Low: " << low_pressure_status << std::endl;
571  return candidates;
572  }
573  else if ( low_pressure_status == 1 && high_pressure_status == 1 )
574  candidates.electron = true;
575 
576  else if ( low_pressure_status == 0 && high_pressure_status == 1 ){
577  candidates.muon = true;
578  candidates.pion = true;
579  }
580 
581  else{ // low, high = 0, 0
582  candidates.proton = true;
583  candidates.kaon = true;
584  }
585  }
586  else if( nominal_momentum == 6. || nominal_momentum == 7. ){
587  if( high_pressure_status == -1 || low_pressure_status == -1 ){
588  std::cout << "At least one Cerenkov status invalid " << std::endl;
589  std::cout << "High: " << high_pressure_status << " Low: " << low_pressure_status << std::endl;
590  return candidates;
591  }
592  else if ( low_pressure_status == 1 && high_pressure_status == 1 ){
593  candidates.electron = true;
594  candidates.muon = true;
595  candidates.pion = true;
596  }
597  else if ( low_pressure_status == 0 && high_pressure_status == 1 )
598  candidates.kaon = true;
599 
600  else // low, high = 0, 0
601  candidates.proton = true;
602  }
603 
604  return candidates;
605 
606 }
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
double protoana::ProtoDUNEBeamlineUtils::GetPosition ( short  theFiber)

Definition at line 242 of file ProtoDUNEBeamlineUtils.cxx.

242  {
243  //Fibers are 1mm in width
244  return 1.*(96 - theFiber) - .5;
245 }
bool protoana::ProtoDUNEBeamlineUtils::HasPerfectBeamMomentum ( art::Event const &  evt) const

Returns true if only 1 fiber is activated in each of the momentum spectrometer monitors

Definition at line 658 of file ProtoDUNEBeamlineUtils.cxx.

658  {
659  return HasPerfectBeamMomentum( GetBeamEvent(evt) );
660 }
protoana::ProtoDUNEBeamlineUtils::HasPerfectBeamMomentum
bool HasPerfectBeamMomentum(art::Event const &evt) const
Definition: ProtoDUNEBeamlineUtils.cxx:658
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
protoana::ProtoDUNEBeamlineUtils::GetBeamEvent
const beam::ProtoDUNEBeamEvent GetBeamEvent(art::Event const &evt) const
Definition: ProtoDUNEBeamlineUtils.cxx:18
bool protoana::ProtoDUNEBeamlineUtils::HasPerfectBeamMomentum ( beam::ProtoDUNEBeamEvent const &  beamEvent) const

Definition at line 662 of file ProtoDUNEBeamlineUtils.cxx.

662  {
663 
664  //Get Active fibers
665  std::vector< short > fibers_p1 = beamEvent.GetActiveFibers( "XBPF022697" );
666  std::vector< short > fibers_p2 = beamEvent.GetActiveFibers( "XBPF022701" );
667  std::vector< short > fibers_p3 = beamEvent.GetActiveFibers( "XBPF022702" );
668 
669  //Get any glitches and remove from the active
670  std::array< short, 192 > glitch_mask = beamEvent.GetFBM( "XBPF022697" ).glitch_mask;
671  int nP1 = 0;
672  for( size_t i = 0; i < fibers_p1.size(); ++i ){
673  if( !glitch_mask[ fibers_p1[i] ] )
674  ++nP1;
675  }
676 
677  glitch_mask = beamEvent.GetFBM( "XBPF022701" ).glitch_mask;
678  int nP2 = 0;
679  for( size_t i = 0; i < fibers_p2.size(); ++i ){
680  if( !glitch_mask[ fibers_p2[i] ] )
681  ++nP2;
682  }
683 
684  glitch_mask = beamEvent.GetFBM( "XBPF022702" ).glitch_mask;
685  int nP3 = 0;
686  for( size_t i = 0; i < fibers_p3.size(); ++i ){
687  if( !glitch_mask[ fibers_p3[i] ] )
688  ++nP3;
689  }
690 
691  return ( nP1 == 1 && nP2 == 1 && nP3 == 1 );
692 }
bool protoana::ProtoDUNEBeamlineUtils::IsGoodBeamlineTrigger ( art::Event const &  evt) const

Returns true if the beamline instrumentation has a good trigger that matches the ProtoDUNE trigger.

Definition at line 649 of file ProtoDUNEBeamlineUtils.cxx.

649  {
650  return IsGoodBeamlineTrigger( GetBeamEvent(evt) );
651 }
protoana::ProtoDUNEBeamlineUtils::IsGoodBeamlineTrigger
bool IsGoodBeamlineTrigger(art::Event const &evt) const
Definition: ProtoDUNEBeamlineUtils.cxx:649
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
protoana::ProtoDUNEBeamlineUtils::GetBeamEvent
const beam::ProtoDUNEBeamEvent GetBeamEvent(art::Event const &evt) const
Definition: ProtoDUNEBeamlineUtils.cxx:18
bool protoana::ProtoDUNEBeamlineUtils::IsGoodBeamlineTrigger ( beam::ProtoDUNEBeamEvent const &  beamEvent) const

Definition at line 653 of file ProtoDUNEBeamlineUtils.cxx.

653  {
654  return (beamEvent.GetTimingTrigger() == 12 && beamEvent.CheckIsMatched());
655 }
std::vector< recob::Track > protoana::ProtoDUNEBeamlineUtils::MakeTracks ( art::Event const &  evt)

Definition at line 72 of file ProtoDUNEBeamlineUtils.cxx.

72  {
73 
74  std::vector< recob::Track > tracks;
75 
76  //Load fibers
77  GetFibers( evt );
78 
79  std::vector<short> HorizUpstreamFibers = ActiveFibers.at(HorizUpstream);
80  std::vector<short> HorizDownstreamFibers = ActiveFibers.at(HorizDownstream);
81  std::vector<short> VertUpstreamFibers = ActiveFibers.at(VertUpstream);
82  std::vector<short> VertDownstreamFibers = ActiveFibers.at(VertDownstream);
83 
84  //Check that the monitors all have at least one hit each
85  bool all_good ;
86  all_good = ( ( HorizUpstreamFibers.size() > 0 ) && ( VertUpstreamFibers.size() > 0 )
87  && ( HorizDownstreamFibers.size() > 0 ) && ( VertDownstreamFibers.size() > 0 ) );
88 
89  if( !all_good ){
90  std::cout << "At least one empty monitor. Producing no track" << std::endl;
91  return tracks;
92  }
93 
94  //Convention: (Horiz, Vert)
95  std::vector< std::pair<short, short> > UpstreamPairedFibers;
96  std::vector< std::pair<short, short> > DownstreamPairedFibers;
97 
98  std::vector< TVector3 > UpstreamPositions;
99  std::vector< TVector3 > DownstreamPositions;
100 
101  for(size_t iH = 0; iH < HorizUpstreamFibers.size(); ++iH){
102 
103  size_t HorizFiber = HorizUpstreamFibers[iH];
104 
105  for(size_t iV = 0; iV < VertUpstreamFibers.size(); ++iV){
106  size_t VertFiber = VertUpstreamFibers[iV];
107 
108  //MF_LOG_DEBUG("BeamEvent") << "Paired: " << HorizFiber << " " << VertFiber << "\n";
109  std::cout << "Paired: " << HorizFiber << " " << VertFiber << std::endl;
110  UpstreamPairedFibers.push_back(std::make_pair(HorizFiber, VertFiber));
111 
112  //If there's 2 adjacent fibers. Skip the next. Could replace this with averaging the position
113  //todo later
114  if (iV < VertUpstreamFibers.size() - 1){
115  if (VertUpstreamFibers[iV] == (VertUpstreamFibers[iV + 1] - 1)) ++iV;
116  }
117  }
118 
119  if (iH < HorizUpstreamFibers.size() - 1){
120  if (HorizUpstreamFibers[iH] == (HorizUpstreamFibers[iH + 1] - 1)) ++iH;
121  }
122  }
123 
124  std::cout << "Upstream " << std::endl;
125 
126  for(size_t i = 0; i < UpstreamPairedFibers.size(); ++i){
127 
128  std::pair<short,short> thePair = UpstreamPairedFibers.at(i);
129 
130  double xPos = GetPosition(thePair.first);
131  double yPos = GetPosition(thePair.second);
132 
133  std::cout << "normal " << xPos << " " << yPos << std::endl;
134  TVector3 posInDet = ConvertMonitorCoordinates(xPos,yPos,0.,fFirstTrackingProfZ);
135  std::cout << posInDet.X() << " " << posInDet.Y() << " " << posInDet.Z() << std::endl;
136  UpstreamPositions.push_back( posInDet );
137  }
138 
139  for(size_t iH = 0; iH < HorizDownstreamFibers.size(); ++iH){
140 
141  size_t HorizFiber = HorizDownstreamFibers[iH];
142 
143  for(size_t iV = 0; iV < VertDownstreamFibers.size(); ++iV){
144  size_t VertFiber = VertDownstreamFibers[iV];
145 
146  //MF_LOG_DEBUG("BeamEvent") << "Paired: " << HorizFiber << " " << VertFiber << "\n";
147  std::cout << "Paired: " << HorizFiber << " " << VertFiber << std::endl;
148  DownstreamPairedFibers.push_back(std::make_pair(HorizFiber, VertFiber));
149 
150  //If there's 2 adjacent fibers. Skip the next. Could replace this with averaging the position
151  //todo later
152  if (iV < VertDownstreamFibers.size() - 1){
153  if (VertDownstreamFibers[iV] == (VertDownstreamFibers[iV + 1] - 1)) ++iV;
154  }
155  }
156 
157  if (iH < HorizDownstreamFibers.size() - 1){
158  if (HorizDownstreamFibers[iH] == (HorizDownstreamFibers[iH + 1] - 1)) ++iH;
159  }
160  }
161 
162  std::cout << "Downstream " << std::endl;
163 
164  for(size_t i = 0; i < DownstreamPairedFibers.size(); ++i){
165 
166  std::pair<short,short> thePair = DownstreamPairedFibers.at(i);
167 
168  double xPos = GetPosition(thePair.first);
169  double yPos = GetPosition(thePair.second);
170 
171  std::cout << "normal " << xPos << " " << yPos << std::endl;
172  TVector3 posInDet = ConvertMonitorCoordinates(xPos,yPos,0.,fSecondTrackingProfZ);
173  std::cout << posInDet.X() << " " << posInDet.Y() << " " << posInDet.Z() << std::endl;
174  DownstreamPositions.push_back( posInDet );
175  }
176 
177  for(size_t iU = 0; iU < UpstreamPositions.size(); ++iU){
178  for(size_t iD = 0; iD < DownstreamPositions.size(); ++iD){
179  std::vector<TVector3> thePoints;
180  thePoints.push_back(UpstreamPositions.at(iU));
181  thePoints.push_back(DownstreamPositions.at(iD));
182 
183  //Now project the last point to the TPC face
184  thePoints.push_back( ProjectToTPC(thePoints[0],thePoints[1]) );
185  std::cout << "Projected: " << thePoints.back().X() << " " << thePoints.back().Y() << " " << thePoints.back().Z() << std::endl;
186 
187 
188  std::vector<TVector3> theMomenta;
189  //Just push back the unit vector for each point
190  theMomenta.push_back( ( DownstreamPositions.at(iD) - UpstreamPositions.at(iU) ).Unit() );
191  theMomenta.push_back( ( DownstreamPositions.at(iD) - UpstreamPositions.at(iU) ).Unit() );
192  theMomenta.push_back( ( DownstreamPositions.at(iD) - UpstreamPositions.at(iU) ).Unit() );
193 
194  recob::Track tempTrack(recob::TrackTrajectory(recob::tracking::convertCollToPoint(thePoints),
195  recob::tracking::convertCollToVector(theMomenta),
196  recob::Track::Flags_t(thePoints.size()), false),
197  0, -1., 0, recob::tracking::SMatrixSym55(), recob::tracking::SMatrixSym55(), 1);
198  tracks.push_back( tempTrack );
199  }
200  }
201 
202  return tracks;
203 }
recob::tracking::SMatrixSym55
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 > > SMatrixSym55
Definition: TrackingTypes.h:85
recob::tracking::convertCollToVector
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:58
protoana::ProtoDUNEBeamlineUtils::ProjectToTPC
TVector3 ProjectToTPC(TVector3, TVector3)
Definition: ProtoDUNEBeamlineUtils.cxx:286
recob::Track::Flags_t
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:68
run_cnn_1class.tracks
tracks
Definition: run_cnn_1class.py:54
protoana::ProtoDUNEBeamlineUtils::GetFibers
void GetFibers(art::Event const &evt)
Definition: ProtoDUNEBeamlineUtils.cxx:42
recob::TrackTrajectory
A trajectory in space reconstructed from hits.
Definition: TrackTrajectory.h:62
recob::tracking::convertCollToPoint
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:55
protoana::ProtoDUNEBeamlineUtils::ConvertMonitorCoordinates
TVector3 ConvertMonitorCoordinates(double, double, double, double)
Definition: ProtoDUNEBeamlineUtils.cxx:247
protoana::ProtoDUNEBeamlineUtils::ActiveFibers
std::map< std::string, std::vector< short > > ActiveFibers
Definition: ProtoDUNEBeamlineUtils.h:227
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
recob::Track
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
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
double protoana::ProtoDUNEBeamlineUtils::MomentumCosTheta ( double  X1,
double  X2,
double  X3 
)

Definition at line 429 of file ProtoDUNEBeamlineUtils.cxx.

429  {
430  double a = (X2*L3 - X3*L2)*cos(fBeamBend)/(L3-L2);
431 
432 
433  double numTerm = (a - X1)*( (L3 - L2)*tan(fBeamBend) + (X3 - X2)*cos(fBeamBend) ) + L1*( L3 - L2 );
434 
435  double denomTerm1, denomTerm2, denom;
436  denomTerm1 = sqrt( L1*L1 + (a - X1)*(a - X1) );
437  denomTerm2 = sqrt( TMath::Power( ( (L3 - L2)*tan(fBeamBend) + (X3 - X2)*cos(fBeamBend) ),2)
438  + TMath::Power( ( (L3 - L2) ),2) );
439  denom = denomTerm1 * denomTerm2;
440 
441  double cosTheta = numTerm/denom;
442  return cosTheta;
443 }
a
const double a
Definition: gUpMuFluxGen.cxx:164
std::vector< double > protoana::ProtoDUNEBeamlineUtils::MomentumSpec ( art::Event const &  evt)

Definition at line 297 of file ProtoDUNEBeamlineUtils.cxx.

297  {
298 
299  GetCurrent( evt );
300 
301  std::cout << "Got current: " << Current << std::endl;
302 
303  std::vector< double > theMomenta;
304 
305  if( Current < .000000001 ){
306  std::cout << "Warning! Low magnet current. Momentum spectrometry invalid. Returning empty momenta vector." << std::endl;
307  return theMomenta;
308  }
309 
310 
311 
312  double LB = mag_P1*fabs(Current);
313  double deltaI = fabs(Current) - mag_P4;
314 
315  if(deltaI>0) LB+= mag_P3*deltaI*deltaI;
316 
317  //Get the active fibers from the upstream tracking XBPF
318  std::vector<short> BProf1Fibers = ActiveFibers.at(BProf1);
319  std::vector<short> BProf2Fibers = ActiveFibers.at(BProf2);
320  std::vector<short> BProf3Fibers = ActiveFibers.at(BProf3);
321 
322 
323  std::cout << BProf1 << " has " << BProf1Fibers.size() << " active fibers" << std::endl;
324  for(size_t i = 0; i < BProf1Fibers.size(); ++i){
325  std::cout << BProf1Fibers[i] << " ";
326  }
327  std::cout << std::endl;
328 
329  std::cout << BProf2 << " has " << BProf2Fibers.size() << " active fibers" << std::endl;
330  for(size_t i = 0; i < BProf2Fibers.size(); ++i){
331  std::cout << BProf2Fibers[i] << " ";
332  }
333  std::cout << std::endl;
334 
335  std::cout << BProf3 << " has " << BProf3Fibers.size() << " active fibers" << std::endl;
336  for(size_t i = 0; i < BProf3Fibers.size(); ++i){
337  std::cout << BProf3Fibers[i] << " ";
338  }
339  std::cout << std::endl;
340 
341 
342  if( (BProf1Fibers.size() < 1) || (BProf2Fibers.size() < 1) || (BProf3Fibers.size() < 1) ){
343  std::cout << "Warning, at least one empty Beam Profiler. Not checking momentum" << std::endl;
344  return theMomenta;
345  }
346 
347 
348  std::cout << "Getting all trio-wise hits" << std::endl;
349  std::cout << "N1,N2,N3 " << BProf1Fibers.size()
350  << " " << BProf2Fibers.size()
351  << " " << BProf3Fibers.size() << std::endl;
352 
353  for(size_t i1 = 0; i1 < BProf1Fibers.size(); ++i1){
354 
355  double x1,x2,x3;
356 
357  x1 = -1.*GetPosition(BProf1Fibers[i1])/1.E3;
358 
359  if (i1 < BProf1Fibers.size() - 1){
360  if (BProf1Fibers[i1] == (BProf1Fibers[i1 + 1] - 1)){
361  //Add .5 mm
362  x1 += .0005;
363  }
364  }
365 
366  for(size_t i2 = 0; i2 < BProf2Fibers.size(); ++i2){
367 
368  x2 = -1.*GetPosition(BProf2Fibers[i2])/1.E3;
369 
370  if (i2 < BProf2Fibers.size() - 1){
371  if (BProf2Fibers[i2] == (BProf2Fibers[i2 + 1] - 1)){
372  //Add .5 mm
373  x2 += .0005;
374  }
375  }
376 
377  for(size_t i3 = 0; i3 < BProf3Fibers.size(); ++i3){
378 
379  std::cout << "\t" << i1 << " " << i2 << " " << i3 << std::endl;
380 
381  x3 = -1.*GetPosition(BProf3Fibers[i3])/1.E3;
382 
383  if (i3 < BProf3Fibers.size() - 1){
384  if (BProf3Fibers[i3] == (BProf3Fibers[i3 + 1] - 1)){
385  //Add .5 mm
386  x3 += .0005;
387  }
388  }
389 
390  //Calibrate the positions of the fibers in the
391  //monitors
392  x1 -= fBProf1Shift*1.e-3;
393  x2 -= fBProf2Shift*1.e-3;
394  x3 -= fBProf3Shift*1.e-3;
395 
396  double cosTheta = MomentumCosTheta(x1,x2,x3);
397  double momentum = 299792458*LB/(1.E9 * acos(cosTheta));
398 
399  theMomenta.push_back(momentum);
400 
401  if (i3 < BProf3Fibers.size() - 1){
402  if (BProf3Fibers[i3] == (BProf3Fibers[i3 + 1] - 1)){
403  //Skip the next
404  ++i3;
405  }
406  }
407 
408  }
409 
410  if (i2 < BProf2Fibers.size() - 1){
411  if (BProf2Fibers[i2] == (BProf2Fibers[i2 + 1] - 1)){
412  //Skip the next
413  ++i2;
414  }
415  }
416  }
417 
418  if (i1 < BProf1Fibers.size() - 1){
419  if (BProf1Fibers[i1] == (BProf1Fibers[i1 + 1] - 1)){
420  //Skip the next
421  ++i1;
422  }
423  }
424  }
425 
426  return theMomenta;
427 }
reco_momentum_tuples.x3
x3
Definition: reco_momentum_tuples.py:84
protoana::ProtoDUNEBeamlineUtils::MomentumCosTheta
double MomentumCosTheta(double, double, double)
Definition: ProtoDUNEBeamlineUtils.cxx:429
reco_momentum_tuples.x2
x2
Definition: reco_momentum_tuples.py:83
reco_momentum_tuples.deltaI
deltaI
Definition: reco_momentum_tuples.py:262
protoana::ProtoDUNEBeamlineUtils::GetCurrent
void GetCurrent(art::Event const &evt)
Definition: ProtoDUNEBeamlineUtils.cxx:27
reco_momentum_tuples.LB
LB
Definition: reco_momentum_tuples.py:261
protoana::ProtoDUNEBeamlineUtils::ActiveFibers
std::map< std::string, std::vector< short > > ActiveFibers
Definition: ProtoDUNEBeamlineUtils.h:227
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
reco_momentum_tuples.momentum
def momentum(x1, x2, x3, scale=1.)
Definition: reco_momentum_tuples.py:198
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
reco_momentum_tuples.x1
x1
Definition: reco_momentum_tuples.py:82
TVector3 protoana::ProtoDUNEBeamlineUtils::ProjectToTPC ( TVector3  firstPoint,
TVector3  secondPoint 
)

Definition at line 286 of file ProtoDUNEBeamlineUtils.cxx.

286  {
287  TVector3 dR = (secondPoint - firstPoint);
288 
289  double deltaZ = -1.*secondPoint.Z();
290  double deltaX = deltaZ * (dR.X() / dR.Z());
291  double deltaY = deltaZ * (dR.Y() / dR.Z());
292 
293  TVector3 lastPoint = secondPoint + TVector3(deltaX, deltaY, deltaZ);
294  return lastPoint;
295 }
void protoana::ProtoDUNEBeamlineUtils::reconfigure ( fhicl::ParameterSet const &  pset)

Definition at line 205 of file ProtoDUNEBeamlineUtils.cxx.

205  {
206  fBeamEventTag = p.get<art::InputTag>("BeamEventTag");
207  fUseCERNCalibSelection = p.get<bool>("UseCERNCalibSelection");
208 
209  //Tracking parameters
210  fBeamX = p.get<double>("BeamX");
211  fBeamY = p.get<double>("BeamY");
212  fBeamZ = p.get<double>("BeamZ");
213 
214  fRotateMonitorXZ = p.get<double>("RotateMonitorXZ");
215  fRotateMonitorYZ = p.get<double>("RotateMonitorYZ");
216  fRotateMonitorYX = p.get<double>("RotateMonitorYX");
217 
218  fFirstTrackingProfZ = p.get<double>("FirstTrackingProfZ");
219  fSecondTrackingProfZ = p.get<double>("SecondTrackingProfZ");
220  fNP04FrontZ = p.get<double>("NP04FrontZ");
221  ////////////////////////
222 
223 
224  //Momentum parameters
225  fBeamBend = p.get<double>("BeamBend");
226  L1 = p.get<double>("L1");
227  L2 = p.get<double>("L2");
228  L3 = p.get<double>("L3");
229 
230  fBProf1Shift = p.get< double >("BProf1Shift");
231  fBProf2Shift = p.get< double >("BProf2Shift");
232  fBProf3Shift = p.get< double >("BProf3Shift");
233  ///////////////////////
234 
235  // Justin's stuff from BeamlineUtils
236  fMomentumScaleFactor = p.get<float>("MomentumScaleFactor");
237  fMomentumOffset = p.get<float>("MomentumOffset"); // GeV/c
238  fTOFScaleFactor = p.get<float>("TOFScaleFactor");
239  fTOFOffset = p.get<float>("TOFOffset"); // ns
240 }
test.p
p
Definition: test.py:223
void protoana::ProtoDUNEBeamlineUtils::RotateMonitorVector ( TVector3 &  vec)

Definition at line 280 of file ProtoDUNEBeamlineUtils.cxx.

280  {
281  vec.RotateX( fRotateMonitorYZ * TMath::Pi()/180. );
282  vec.RotateY( fRotateMonitorXZ * TMath::Pi()/180. );
283 }

Member Data Documentation

std::map< std::string, std::vector< short > > protoana::ProtoDUNEBeamlineUtils::ActiveFibers
private

Definition at line 227 of file ProtoDUNEBeamlineUtils.h.

std::vector< std::string > protoana::ProtoDUNEBeamlineUtils::AllDevices
private
Initial value:
= {
"XBPF022707",
"XBPF022708",
"XBPF022716",
"XBPF022717",
"XBPF022697",
"XBPF022701",
"XBPF022702"
}

Definition at line 235 of file ProtoDUNEBeamlineUtils.h.

std::string protoana::ProtoDUNEBeamlineUtils::BProf1 = "XBPF022697"
private

Definition at line 254 of file ProtoDUNEBeamlineUtils.h.

std::string protoana::ProtoDUNEBeamlineUtils::BProf2 = "XBPF022701"
private

Definition at line 255 of file ProtoDUNEBeamlineUtils.h.

std::string protoana::ProtoDUNEBeamlineUtils::BProf3 = "XBPF022702"
private

Definition at line 256 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::c = 299792458.
private

Definition at line 295 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::Current
private

Definition at line 230 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fBeamBend
private

Definition at line 276 of file ProtoDUNEBeamlineUtils.h.

art::InputTag protoana::ProtoDUNEBeamlineUtils::fBeamEventTag
private

Definition at line 222 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fBeamX
private

Definition at line 273 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fBeamY
private

Definition at line 273 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fBeamZ
private

Definition at line 273 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fBProf1Shift
private

Definition at line 278 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fBProf2Shift
private

Definition at line 278 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fBProf3Shift
private

Definition at line 278 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fFirstTrackingProfZ
private

Definition at line 270 of file ProtoDUNEBeamlineUtils.h.

float protoana::ProtoDUNEBeamlineUtils::fMomentumOffset
private

Definition at line 302 of file ProtoDUNEBeamlineUtils.h.

float protoana::ProtoDUNEBeamlineUtils::fMomentumScaleFactor
private

Definition at line 301 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fNP04FrontZ
private

Definition at line 272 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fRotateMonitorXZ
private

Definition at line 266 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fRotateMonitorYX
private

Definition at line 268 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fRotateMonitorYZ
private

Definition at line 267 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::fSecondTrackingProfZ
private

Definition at line 271 of file ProtoDUNEBeamlineUtils.h.

const double protoana::ProtoDUNEBeamlineUtils::fTOFDist = 28.575
private

Definition at line 296 of file ProtoDUNEBeamlineUtils.h.

float protoana::ProtoDUNEBeamlineUtils::fTOFOffset
private

Definition at line 304 of file ProtoDUNEBeamlineUtils.h.

float protoana::ProtoDUNEBeamlineUtils::fTOFScaleFactor
private

Definition at line 303 of file ProtoDUNEBeamlineUtils.h.

bool protoana::ProtoDUNEBeamlineUtils::fUseCERNCalibSelection
private

Definition at line 224 of file ProtoDUNEBeamlineUtils.h.

std::string protoana::ProtoDUNEBeamlineUtils::HorizDownstream = "XBPF022716"
private

Definition at line 250 of file ProtoDUNEBeamlineUtils.h.

std::string protoana::ProtoDUNEBeamlineUtils::HorizUpstream = "XBPF022707"
private

Definition at line 248 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::L1
private

Definition at line 277 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::L2
private

Definition at line 277 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::L3
private

Definition at line 277 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::mag_P1 = 5.82044830e-3
private

Definition at line 281 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::mag_P3 = -4.68880000e-6
private

Definition at line 283 of file ProtoDUNEBeamlineUtils.h.

double protoana::ProtoDUNEBeamlineUtils::mag_P4 = 324.573967
private

Definition at line 284 of file ProtoDUNEBeamlineUtils.h.

TVector3 protoana::ProtoDUNEBeamlineUtils::MonitorBasisX
private

Definition at line 260 of file ProtoDUNEBeamlineUtils.h.

TVector3 protoana::ProtoDUNEBeamlineUtils::MonitorBasisY
private

Definition at line 261 of file ProtoDUNEBeamlineUtils.h.

TVector3 protoana::ProtoDUNEBeamlineUtils::MonitorBasisZ
private

Definition at line 262 of file ProtoDUNEBeamlineUtils.h.

std::map< int, double > protoana::ProtoDUNEBeamlineUtils::particle_mass
private
Initial value:
= {
{kElectron, .0005109989},
{kMuon, .1056583745},
{kPion, .13957018 },
{kKaon, .493677 },
{kProton, .9382720813},
{kDeuteron, 2.013553213}
}

Definition at line 286 of file ProtoDUNEBeamlineUtils.h.

bool protoana::ProtoDUNEBeamlineUtils::rotated = false
private

Definition at line 264 of file ProtoDUNEBeamlineUtils.h.

std::string protoana::ProtoDUNEBeamlineUtils::VertDownstream = "XBPF022717"
private

Definition at line 251 of file ProtoDUNEBeamlineUtils.h.

std::string protoana::ProtoDUNEBeamlineUtils::VertUpstream = "XBPF022708"
private

Definition at line 249 of file ProtoDUNEBeamlineUtils.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11