ProtoDUNEBeam_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       ProtoDUNEBeam
// Module Type: producer
// File:        ProtoDUNEBeam_module.cc
//
// Generated at Thu Nov 17 11:20:31 2016 by Leigh Howard Whitehead,42 3-039,+41227672470, using artmod
// from cetpkgsupport v1_11_00.
////////////////////////////////////////////////////////////////////////
// Modified by Caroline Zhang with inputs from Karl Wharburton  
// August 2017 for beam simulation storage
// Email: carolineligezhang@gmail.com
////////////////////////////////////////////////////////////////////////
// Modified by Pablo and Leigh H. Whitehead  
// July 2018 for redoing beam simulation storage and definition of 
// absent Cherenkov detectors
// Email: pablo.fer@cern.ch
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSetRegistry.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "cetlib_except/exception.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "larcore/Geometry/Geometry.h"
#include "larcoreobj/SummaryData/RunData.h"
#include <memory>
#include <string>
#include <map>
#include <utility>
#include <vector>
// Added for ProtoDUNE beam simulation storage
//#include "lardataobj/Simulation/ProtoDUNEbeamsim.h"
#include "dunesim/EventGenerator/ProtoDUNEbeamDataProducts/ProtoDUNEbeamsim.h"
#include "dunecore/DuneObj/ProtoDUNEBeamEvent.h"
#include "dunecore/DuneObj/ProtoDUNEBeamSpill.h"
#include "lardataobj/RecoBase/TrackingTypes.h"
#include "lardataobj/RecoBase/TrackTrajectory.h"
#include "lardataobj/RecoBase/Track.h"
#include "dunesim/EventGenerator/ProtoDUNEbeamDataProducts/ProtoDUNEBeamInstrument.h"
//#include "dunesim/EventGenerator/ProtoDUNEbeamTPCmatching/ProtoDUNEbeammatch.h"
//#include "dunesim/EventGenerator/ProtoDUNEbeamTPCmatching/ProtoDUNEBeamToF.h"
#include "lardata/Utilities/AssociationUtil.h"
// art extensions
#include "nurandom/RandomUtils/NuRandomService.h"
#include "art_root_io/TFileService.h"
#include <TFile.h>
#include <TTree.h>
#include <TVector3.h>
#include <TLorentzVector.h>
#include <TDatabasePDG.h>
#include <TParticlePDG.h>
#include "TSystem.h"
#include "CLHEP/Random/RandFlat.h"
#include "ifdh.h"
#include <sys/stat.h>

namespace evgen{

    class ProtoDUNEBeam;
 
    class ProtoDUNEBeam : public art::EDProducer {
    public:
        explicit ProtoDUNEBeam(fhicl::ParameterSet const & p);
        // The destructor generated by the compiler is fine for classes
        // without bare pointers or other resource use.
        ~ProtoDUNEBeam();
        
        // Plugins should not be copied or assigned.
        ProtoDUNEBeam(ProtoDUNEBeam const &) = delete;
        ProtoDUNEBeam(ProtoDUNEBeam &&) = delete;
        ProtoDUNEBeam & operator = (ProtoDUNEBeam const &) = delete;
        ProtoDUNEBeam & operator = (ProtoDUNEBeam &&) = delete;
        
        // Required functions.
        void produce(art::Event & e) override;
        void beginJob() override;
        void beginRun(art::Run& run) override;
        void endJob() override;


        // Convenience struct to encapsulate each spill
        // Contains the good particle and all backgrounds
        struct ProtoFullSpill {

          ProtoFullSpill(int event, int track, float time, int id){
            fGoodEvent = event;
            fGoodTrack = track;
            fGoodTime = time;
            fGoodIndex = id;
          };

          // Good particle information
          int fGoodEvent; // Beam simulation event number
          int fGoodTrack; // Beam simulation track number
          float fGoodTime;  // Time of beam event
          int fGoodIndex; // Index of the good particle in the good particle tree

          // All of the beam events and tracks
          std::map<int,std::vector<int> > fAllSpillTracks; 
        };
        
    private:
        
        std::vector<ProtoFullSpill> fAllSpills;

        // A list of good events and an index for it.
        unsigned int fCurrentGoodEvent;
        std::vector<int> fGoodEventList;
        
        // Calculate how many overlay events we need.
        void CalculateNOverlays();
        
        // Check if a given beam event is close enough to a good particle event to be useful.
        int IsOverlayEvent(int event, int nOverlay);
        std::vector<int> GetAllOverlays(int event, int nOverlay);
        
        // Fill the above maps and vector.
        void FillParticleMaps();
        
        // Generate a true event based on a single entry from the input tree.
        void GenerateTrueEvent(simb::MCTruth &mcTruth, std::vector<sim::ProtoDUNEbeamsim> &beamsimcol, beam::ProtoDUNEBeamEvent & beamEvent);
        
        // Handle root files from beam instrumentation group
        void OpenInputFile();
        
        // Generate a TLorentzVector for position making sure we get the
        // coordinates as we need them.
        TLorentzVector ConvertCoordinates(float x, float y, float z, float t);
        
        // Make the momentum vector, rotating as required.
        TLorentzVector MakeMomentumVector(float px, float py, float pz, int pdg, bool shifts);
        TLorentzVector MakeMomentumVector(const TVector3 &mom, int pdg, bool shifts);

        // We need to rotate the beam monitor coordinates into the detector frame
        TLorentzVector ConvertBeamMonitorCoordinates(float x, float y, float z, float t, float offset);

        //New for making tracks
        TVector3 ConvertProfCoordinates(double x, double y, double z, double zOffset);
        TVector3 ProjectToTPC(TVector3 firstPoint, TVector3 secondPoint);
        double GetPosition( short fiber );
        void MakeTracks( beam::ProtoDUNEBeamEvent & beamEvent );
        void MomentumSpectrometer( beam::ProtoDUNEBeamEvent & beamEvent );
        double MomentumCosTheta( double, double, double );


        TVector3 ConvertBeamMonitorMomentumVec(float px, float py, float pz);
        // Setup the beam monitor basis vectors in detector coordinates
        void BeamMonitorBasisVectors();
        // Apply the rotation
        void RotateMonitorVector(TVector3 &vec);

        void SetBeamEvent(beam::ProtoDUNEBeamEvent & beamevt);
        beam::FBM MakeFiberMonitor( float pos );

 
        // Background particles need to be fired from an upstream position
        TVector3 GetBackgroundPosition(float x, float y, float z, float px, float py, float pz);
       
        CLHEP::RandFlat fFlatRnd;

        std::string fFileName;
        std::string fGoodParticleTreeName;
        std::string fAllParticlesTreeName;
        
        // The current event number. Ideally this could be an unsigned int,
        // but we will need to compare it to some ints later on.
        int fEventNumber;
        
        // Let the user define the event to start at
        int fStartEvent;
        
        TFile* fInputFile;
        // Input file provides a TTree that we need to read.
        TTree* fGoodParticleTree;
        TTree* fAllParticlesTree;

        // To make sure we can fire the GoodParticles from the correct place        
        Float_t fGoodNP04front_x;
        Float_t fGoodNP04front_y;
        Float_t fGoodNP04front_z;
        Float_t fGoodNP04front_t;
        Float_t fGoodNP04front_Px;
        Float_t fGoodNP04front_Py;
        Float_t fGoodNP04front_Pz;
        Float_t fGoodNP04front_PDGid;
        Float_t fGoodNP04front_EventID;
        Float_t fGoodNP04front_TrackID;

        // For the TOF part
        Float_t fGoodTOF1_x;
        Float_t fGoodTOF1_y;
        Float_t fGoodTOF1_z;
        Float_t fGoodTOF1_t;
        Float_t fGoodTOF1_Px;
        Float_t fGoodTOF1_Py;
        Float_t fGoodTOF1_Pz;
        Float_t fGoodTOF1_PDGid;
        Float_t fGoodTOF1_EventID;
        Float_t fGoodTOF1_TrackID;

        //For the TRIG part
        Float_t fGoodTRIG1_x;
        Float_t fGoodTRIG1_y;
        Float_t fGoodTRIG1_z;
        Float_t fGoodTRIG1_t;
        Float_t fGoodTRIG1_Px;
        Float_t fGoodTRIG1_Py;
        Float_t fGoodTRIG1_Pz;
        Float_t fGoodTRIG1_PDGid;
        Float_t fGoodTRIG1_EventID;
        Float_t fGoodTRIG1_TrackID;

        Float_t fGoodTRIG2_x;
        Float_t fGoodTRIG2_y;
        Float_t fGoodTRIG2_z;
        Float_t fGoodTRIG2_t;
        Float_t fGoodTRIG2_Px;
        Float_t fGoodTRIG2_Py;
        Float_t fGoodTRIG2_Pz;
        Float_t fGoodTRIG2_PDGid;
        Float_t fGoodTRIG2_EventID;
        Float_t fGoodTRIG2_TrackID;

        //For the BPROF part
        Float_t fGoodBPROF1_x;
        Float_t fGoodBPROF1_y;
        Float_t fGoodBPROF1_z;
        Float_t fGoodBPROF1_t;
        Float_t fGoodBPROF1_Px;
        Float_t fGoodBPROF1_Py;
        Float_t fGoodBPROF1_Pz;
        Float_t fGoodBPROF1_PDGid;
        Float_t fGoodBPROF1_EventID;
        Float_t fGoodBPROF1_TrackID;

        Float_t fGoodBPROF2_x;
        Float_t fGoodBPROF2_y;
        Float_t fGoodBPROF2_z;
        Float_t fGoodBPROF2_t;
        Float_t fGoodBPROF2_Px;
        Float_t fGoodBPROF2_Py;
        Float_t fGoodBPROF2_Pz;
        Float_t fGoodBPROF2_PDGid;
        Float_t fGoodBPROF2_EventID;
        Float_t fGoodBPROF2_TrackID;

        Float_t fGoodBPROF3_x;
        Float_t fGoodBPROF3_y;
        Float_t fGoodBPROF3_z;
        Float_t fGoodBPROF3_t;
        Float_t fGoodBPROF3_Px;
        Float_t fGoodBPROF3_Py;
        Float_t fGoodBPROF3_Pz;
        Float_t fGoodBPROF3_PDGid;
        Float_t fGoodBPROF3_EventID;
        Float_t fGoodBPROF3_TrackID;

        Float_t fGoodBPROF4_x;
        Float_t fGoodBPROF4_y;
        Float_t fGoodBPROF4_z;
        Float_t fGoodBPROF4_t;
        Float_t fGoodBPROF4_Px;
        Float_t fGoodBPROF4_Py;
        Float_t fGoodBPROF4_Pz;
        Float_t fGoodBPROF4_PDGid;
        Float_t fGoodBPROF4_EventID;
        Float_t fGoodBPROF4_TrackID;

        Float_t fGoodBPROFEXT_x;
        Float_t fGoodBPROFEXT_y;
        Float_t fGoodBPROFEXT_z;
        Float_t fGoodBPROFEXT_t;
        Float_t fGoodBPROFEXT_Px;
        Float_t fGoodBPROFEXT_Py;
        Float_t fGoodBPROFEXT_Pz;
        Float_t fGoodBPROFEXT_PDGid;
        Float_t fGoodBPROFEXT_EventID;
        Float_t fGoodBPROFEXT_TrackID;

        TTree * fRecoTree;
        float   fXBPF697_x;
        float   fXBPF701_x;
        float   fXBPF702_x;

        float   fXBPF697_rx;
        float   fXBPF701_rx;
        float   fXBPF702_rx;

        float   fXBPF697_p;
        float   fXBPF701_p;
        float   fXBPF702_p;

        float   fXBPF697_f;
        float   fXBPF701_f;
        float   fXBPF702_f;
        float   fReco_p;

        float   fReco_tof;

        int     fNP04_PDG;
        float   fNP04front_p;

        float   fXBPF716_x;
        float   fXBPF717_y;
        float   fXBPF707_x;
        float   fXBPF708_y;

        float   fXBPF716_rx;
        float   fXBPF717_ry;
        float   fXBPF707_rx;
        float   fXBPF708_ry;

        float   fXBPF716_f;
        float   fXBPF717_f;
        float   fXBPF707_f;
        float   fXBPF708_f;

        float   fTrueFront_x;
        float   fTrueFront_y;
        float   fTrueFront_z;

        float   fTrueFront_Px;
        float   fTrueFront_Py;
        float   fTrueFront_Pz;

        float   fRecoFront_x;
        float   fRecoFront_y;
        float   fRecoFront_z;

        // Members we need to extract from the tree
        float fX, fY, fZ;
        float fPx, fPy, fPz;
        float fPDG; // Input tree has all floats

        int fTrueID;

        // Event and TrackID for good particle tree
//        float fBeamEvent;
//        float fTrackID;

        // Same for all particles
        float fAllEventID;
        float fAllTrackID;
        float fAllParentID;

        // We need two times: the trigger time, and the time at the entry point
        // to the TPC where we generate the event.
        float fEntryT;
        
        // Define the coordinate transform from the beam frame to the detector frame
        float fBeamX;
        float fBeamY;
        float fBeamZ;
        float fBeamThetaShift;
        float fBeamPhiShift;
        float fRotateXZ;
        float fRotateYZ;
        // Rotate the beam monitor coordinate system (those after the last bending magnet)
        float fRotateMonitorXZ;
        float fRotateMonitorYZ;
        // The three beam monitor basis vectors in the detector coordinate system
        TVector3 fBMBasisX; 
        TVector3 fBMBasisY; 
        TVector3 fBMBasisZ; 
        // The z positions of the important elements along the beam direction
        float fBPROFEXTPos;
        float fBPROF4Pos;
        float fNP04frontPos;
       
        // Parameters from the .fcl file to deal with overlaying events
        float fIntensity; // Number of interactions on the secondary target per SPS spill
        float fReadoutWindow; // Readout window (needs to match the values used in the simulation) in milliseconds
        float fBeamSpillLength; // The SPS spill length in seconds

// Beam monitors resolutions
        float fT_Resolution;
        float fPos_Resolution;
        float fCh_Efficiency;

        float fLB;
/*        
        float fMagP1;
        float fMagP3;
        float fMagP4;
        float fCurrent;
*/        
        float fL1;
        float fL2;
        float fL3;
        float fBeamBend;

        float fLMag;
        float fNominalP;
        float fB;

        bool fSaveRecoTree;

        
        // Number of beam interactions to overlay.
        int fOverlays;
        
        ifdh_ns::ifdh* fIFDH;
    };
}

// Create the random number generator
namespace {
  std::string const instanceName = "protoDUNEBeam";
}


//---------------------------------------------------------------------------------
//----------------------------------------constructors-----------------------------
evgen::ProtoDUNEBeam::ProtoDUNEBeam(fhicl::ParameterSet const & pset)
  : EDProducer{pset}
    // now create the engine (for example, use art); seed will be set
    // by calling declareEngine
  , fFlatRnd(createEngine(art::ServiceHandle<rndm::NuRandomService>{}->declareEngine(instanceName),
                          "HepJamesRandom", instanceName))
{
    // Call appropriate produces<>() functions here.
    produces< std::vector<simb::MCTruth> >();
    produces<std::vector<sim::ProtoDUNEbeamsim>>();
    produces< sumdata::RunData, art::InRun >();
    produces< std::vector< beam::ProtoDUNEBeamEvent > >();
    //produces< art::Assns<sim::ProtoDUNEbeamsim, simb::MCTruth>>();
    // File reading variable initialisations
    fFileName = pset.get< std::string>("FileName");
    fGoodParticleTreeName = pset.get< std::string>("GoodParticleTreeName");
    fAllParticlesTreeName = pset.get< std::string>("AllParticlesTreeName");
    std::cout << "All particles tree name = " << fAllParticlesTreeName << std::endl;
    // Intensity variables
    fIntensity = pset.get<float>("Intensity");
    fReadoutWindow = pset.get<float>("ReadoutWindow");
    fBeamSpillLength = pset.get<float>("BeamSpillLength");

// Beam monitors resolutions
    fT_Resolution = pset.get<float>("T_Resolution");
    fPos_Resolution = pset.get<float>("Pos_Resolution");
    fCh_Efficiency = pset.get<float>("Ch_Efficiency");
    
    // See if the user wants to start at an event other than zero.
    fStartEvent = pset.get<int>("StartEvent");
    
    // Or maybe there was --nskip specified in the command line or skipEvents in FHiCL?
    for (auto const & p : fhicl::ParameterSetRegistry::get())
    {
        auto const& ps = p.second;
        if (ps.has_key("source") && ps.has_key("source.skipEvents"))
        {
            fStartEvent += ps.get<int>("source.skipEvents");
            break; // take the first occurence
        } // no "else", if parameter not found, then just don't change anything
    }
    // ...and if there is -e option or firstEvent in FHiCL, this add up to the no. of events to skip.
    for (auto const & p : fhicl::ParameterSetRegistry::get())
    {
        auto const& ps = p.second;
        if (ps.has_key("source") && ps.has_key("source.firstEvent"))
        {
            int fe = ps.get<int>("source.firstEvent") - 1; // events base index is 1
            if (fe > 0) fStartEvent += fe;
            break; // take the first occurence
        } // no "else", if parameter not found, then just don't change anything
    }
    mf::LogInfo("ProtoDUNEBeam") << "Skip " << fStartEvent << " first events from the input file.";
    
    fEventNumber = 0;
    
    // Coordinate transform
    fBeamX = pset.get<float>("BeamX");
    fBeamY = pset.get<float>("BeamY");
    fBeamZ = pset.get<float>("BeamZ");
    fBeamThetaShift = pset.get<float>("BeamThetaShift",0.0);
    fBeamPhiShift   = pset.get<float>("BeamPhiShift",0.0);
    fRotateXZ = pset.get<float>("RotateXZ"); // Only use rotation if the above aren't defined
    fRotateYZ = pset.get<float>("RotateYZ");
    
    fRotateMonitorXZ = pset.get<float>("RotateMonitorXZ");
    fRotateMonitorYZ = pset.get<float>("RotateMonitorYZ");
    fBPROFEXTPos     = pset.get<float>("BPROFEXTPosZ");
    fBPROF4Pos       = pset.get<float>("BPROF4PosZ");
    fNP04frontPos    = pset.get<float>("NP04frontPosZ");
    // Setup the beam monitor basis vectors
    BeamMonitorBasisVectors();   

    // Initialise the input file and tree to be null.
    fInputFile = 0x0;
    fGoodParticleTree = 0x0;
    fAllParticlesTree = 0x0;
    fIFDH = 0;
    
    fCurrentGoodEvent = 0;

    // For momentum spectrometer
/*
    fCurrent = pset.get<float>("Current");
    fMagP1 = pset.get<float>("MagP1");
    fMagP3 = pset.get<float>("MagP3");
    fMagP4 = pset.get<float>("MagP4");
*/    
    fL1 = pset.get<float>("L1");
    fL2 = pset.get<float>("L2");
    fL3 = pset.get<float>("L3");
    fBeamBend = pset.get<float>("BeamBend");

    //New values for momentum spectrometer
    fLMag = pset.get<float>("LMag");
    fB    = pset.get<float>("B");
    fNominalP = pset.get<float>("NominalP");


/*
    fLB = fMagP1*fabs(fCurrent);
    float deltaI = fabs(fCurrent) - fMagP4;
    if(deltaI>0) fLB += fMagP3*deltaI*deltaI;

    std::cout << "Old LB: " << fLB << std::endl;
*/
    fLB = fB * fLMag * fNominalP / 7.;
//    std::cout << "New LB: " << fLB << std::endl;
    
    fSaveRecoTree = pset.get<bool>("SaveRecoTree");

    
    // Make sure we use ifdh to open the beam input file.
    OpenInputFile();
}



//-----------------------------default destructor----------------------------------
//------------------------------------------------------------------------------------
evgen::ProtoDUNEBeam::~ProtoDUNEBeam()
{
    fIFDH->cleanup();
    
}

//-------------------------------------------------------------------------------------
void evgen::ProtoDUNEBeam::beginJob(){

    art::ServiceHandle<art::TFileService> tfs;
    
    fInputFile = new TFile(fFileName.c_str(),"READ");
    // Check we have the file
    if(fInputFile == 0x0){
        throw cet::exception("ProtoDUNEBeam") << "Input file " << fFileName << " cannot be read.\n";
    }
    
    fGoodParticleTree = (TTree*)fInputFile->Get(fGoodParticleTreeName.c_str());
    // Check we have the tree
    if(fGoodParticleTree == 0x0){
        throw cet::exception("ProtoDUNEBeam") << "Input tree " << fGoodParticleTreeName << " cannot be read.\n";
    }
    
    fAllParticlesTree = (TTree*)fInputFile->Get(fAllParticlesTreeName.c_str());
    // Check we have the tree
    if(fAllParticlesTree == 0x0){
        throw cet::exception("ProtoDUNEBeam") << "Input tree " << fAllParticlesTreeName << " cannot be read.\n";
    }
    std::cout << "All particle tree " << fAllParticlesTreeName << " has " << fAllParticlesTree->GetEntries() << " entries" << std::endl;
   
    // We need different bits of information for different particles.
    // The GoodParticles should be fired from the detector front (NP04front)
    // Background particles should be fired from the first monitor after the final bending magnet in order
    // to ensure those that should hit the CRT will do so (BPROFEXT)
 
    // Since this is technically an ntuple, all objects are floats
    // Position four-vector components
    fAllParticlesTree->SetBranchAddress("x",&fX);
    fAllParticlesTree->SetBranchAddress("y",&fY);
    fAllParticlesTree->SetBranchAddress("z",&fZ);
    fAllParticlesTree->SetBranchAddress("t",&fEntryT);
    // Momentum components
    fAllParticlesTree->SetBranchAddress("Px",&fPx);
    fAllParticlesTree->SetBranchAddress("Py",&fPy);
    fAllParticlesTree->SetBranchAddress("Pz",&fPz);
    // PDG code
    fAllParticlesTree->SetBranchAddress("PDGid",&fPDG);
    // Event and track number
    fAllParticlesTree->SetBranchAddress("EventID",&fAllEventID);
    fAllParticlesTree->SetBranchAddress("TrackID",&fAllTrackID);
    fAllParticlesTree->SetBranchAddress("ParentID", &fAllParentID);
    
    // We only need the trigger time and event number from the good particle tree.
    // The good particle tree variable should match the names of the other trees
    std::string namePrefix = fAllParticlesTreeName.substr(fAllParticlesTreeName.find_last_of("\\/")+1,std::string::npos);

    std::cout << "Name prefix for good particle tree = " << namePrefix << std::endl;

//    fGoodParticleTree->SetBranchAddress((namePrefix+"_EventID").c_str(),&fBeamEvent);
//    fGoodParticleTree->SetBranchAddress((namePrefix+"_TrackID").c_str(),&fTrackID);
 
    ////////************added by Caroline for beam simulation storage for good particles ***************//////////////
    fGoodParticleTree->SetBranchAddress("NP04front_x",&fGoodNP04front_x);
    fGoodParticleTree->SetBranchAddress("NP04front_y",&fGoodNP04front_y);
    fGoodParticleTree->SetBranchAddress("NP04front_z",&fGoodNP04front_z);
    fGoodParticleTree->SetBranchAddress("NP04front_t",&fGoodNP04front_t);
    fGoodParticleTree->SetBranchAddress("NP04front_Px",&fGoodNP04front_Px);
    fGoodParticleTree->SetBranchAddress("NP04front_Py",&fGoodNP04front_Py);
    fGoodParticleTree->SetBranchAddress("NP04front_Pz",&fGoodNP04front_Pz);
    fGoodParticleTree->SetBranchAddress("NP04front_PDGid",&fGoodNP04front_PDGid);
    fGoodParticleTree->SetBranchAddress("NP04front_EventID",&fGoodNP04front_EventID);
    fGoodParticleTree->SetBranchAddress("NP04front_TrackID",&fGoodNP04front_TrackID);
    
    // add more lines for the TOF part
    fGoodParticleTree->SetBranchAddress("TOF1_x",&fGoodTOF1_x);
    fGoodParticleTree->SetBranchAddress("TOF1_y",&fGoodTOF1_y);
    fGoodParticleTree->SetBranchAddress("TOF1_z",&fGoodTOF1_z);
    fGoodParticleTree->SetBranchAddress("TOF1_t",&fGoodTOF1_t);
    fGoodParticleTree->SetBranchAddress("TOF1_Px",&fGoodTOF1_Px);
    fGoodParticleTree->SetBranchAddress("TOF1_Py",&fGoodTOF1_Py);
    fGoodParticleTree->SetBranchAddress("TOF1_Pz",&fGoodTOF1_Pz);
    fGoodParticleTree->SetBranchAddress("TOF1_PDGid",&fGoodTOF1_PDGid);
    fGoodParticleTree->SetBranchAddress("TOF1_EventID",&fGoodTOF1_EventID);
    fGoodParticleTree->SetBranchAddress("TOF1_TrackID",&fGoodTOF1_TrackID);

    // add more lines for the TRIG part
    fGoodParticleTree->SetBranchAddress("TRIG1_x",&fGoodTRIG2_x);
    fGoodParticleTree->SetBranchAddress("TRIG1_y",&fGoodTRIG2_y);
    fGoodParticleTree->SetBranchAddress("TRIG1_z",&fGoodTRIG2_z);
    fGoodParticleTree->SetBranchAddress("TRIG1_t",&fGoodTRIG2_t);
    fGoodParticleTree->SetBranchAddress("TRIG1_Px",&fGoodTRIG2_Px);
    fGoodParticleTree->SetBranchAddress("TRIG1_Py",&fGoodTRIG2_Py);
    fGoodParticleTree->SetBranchAddress("TRIG1_Pz",&fGoodTRIG2_Pz);
    fGoodParticleTree->SetBranchAddress("TRIG1_PDGid",&fGoodTRIG2_PDGid);
    fGoodParticleTree->SetBranchAddress("TRIG1_EventID",&fGoodTRIG2_EventID);
    fGoodParticleTree->SetBranchAddress("TRIG1_TrackID",&fGoodTRIG2_TrackID);
    
    fGoodParticleTree->SetBranchAddress("TRIG2_x",&fGoodTRIG2_x);
    fGoodParticleTree->SetBranchAddress("TRIG2_y",&fGoodTRIG2_y);
    fGoodParticleTree->SetBranchAddress("TRIG2_z",&fGoodTRIG2_z);
    fGoodParticleTree->SetBranchAddress("TRIG2_t",&fGoodTRIG2_t);
    fGoodParticleTree->SetBranchAddress("TRIG2_Px",&fGoodTRIG2_Px);
    fGoodParticleTree->SetBranchAddress("TRIG2_Py",&fGoodTRIG2_Py);
    fGoodParticleTree->SetBranchAddress("TRIG2_Pz",&fGoodTRIG2_Pz);
    fGoodParticleTree->SetBranchAddress("TRIG2_PDGid",&fGoodTRIG2_PDGid);
    fGoodParticleTree->SetBranchAddress("TRIG2_EventID",&fGoodTRIG2_EventID);
    fGoodParticleTree->SetBranchAddress("TRIG2_TrackID",&fGoodTRIG2_TrackID);

    //add more lines for the BPROF part
    fGoodParticleTree->SetBranchAddress("BPROF1_x",&fGoodBPROF1_x);
    fGoodParticleTree->SetBranchAddress("BPROF1_y",&fGoodBPROF1_y);
    fGoodParticleTree->SetBranchAddress("BPROF1_z",&fGoodBPROF1_z);
    fGoodParticleTree->SetBranchAddress("BPROF1_t",&fGoodBPROF1_t);
    fGoodParticleTree->SetBranchAddress("BPROF1_Px",&fGoodBPROF1_Px);
    fGoodParticleTree->SetBranchAddress("BPROF1_Py",&fGoodBPROF1_Py);
    fGoodParticleTree->SetBranchAddress("BPROF1_Pz",&fGoodBPROF1_Pz);
    fGoodParticleTree->SetBranchAddress("BPROF1_PDGid",&fGoodBPROF1_PDGid);
    fGoodParticleTree->SetBranchAddress("BPROF1_EventID",&fGoodBPROF1_EventID);
    fGoodParticleTree->SetBranchAddress("BPROF1_TrackID",&fGoodBPROF1_TrackID);

    fGoodParticleTree->SetBranchAddress("BPROF2_x",&fGoodBPROF2_x);
    fGoodParticleTree->SetBranchAddress("BPROF2_y",&fGoodBPROF2_y);
    fGoodParticleTree->SetBranchAddress("BPROF2_z",&fGoodBPROF2_z);
    fGoodParticleTree->SetBranchAddress("BPROF2_t",&fGoodBPROF2_t);
    fGoodParticleTree->SetBranchAddress("BPROF2_Px",&fGoodBPROF2_Px);
    fGoodParticleTree->SetBranchAddress("BPROF2_Py",&fGoodBPROF2_Py);
    fGoodParticleTree->SetBranchAddress("BPROF2_Pz",&fGoodBPROF2_Pz);
    fGoodParticleTree->SetBranchAddress("BPROF2_PDGid",&fGoodBPROF2_PDGid);
    fGoodParticleTree->SetBranchAddress("BPROF2_EventID",&fGoodBPROF2_EventID);
    fGoodParticleTree->SetBranchAddress("BPROF2_TrackID",&fGoodBPROF2_TrackID);

    fGoodParticleTree->SetBranchAddress("BPROF3_x",&fGoodBPROF3_x);
    fGoodParticleTree->SetBranchAddress("BPROF3_y",&fGoodBPROF3_y);
    fGoodParticleTree->SetBranchAddress("BPROF3_z",&fGoodBPROF3_z);
    fGoodParticleTree->SetBranchAddress("BPROF3_t",&fGoodBPROF3_t);
    fGoodParticleTree->SetBranchAddress("BPROF3_Px",&fGoodBPROF3_Px);
    fGoodParticleTree->SetBranchAddress("BPROF3_Py",&fGoodBPROF3_Py);
    fGoodParticleTree->SetBranchAddress("BPROF3_Pz",&fGoodBPROF3_Pz);
    fGoodParticleTree->SetBranchAddress("BPROF3_PDGid",&fGoodBPROF3_PDGid);
    fGoodParticleTree->SetBranchAddress("BPROF3_EventID",&fGoodBPROF3_EventID);
    fGoodParticleTree->SetBranchAddress("BPROF3_TrackID",&fGoodBPROF3_TrackID);

    fGoodParticleTree->SetBranchAddress("BPROF4_x",&fGoodBPROF4_x);
    fGoodParticleTree->SetBranchAddress("BPROF4_y",&fGoodBPROF4_y);
    fGoodParticleTree->SetBranchAddress("BPROF4_z",&fGoodBPROF4_z);
    fGoodParticleTree->SetBranchAddress("BPROF4_t",&fGoodBPROF4_t);
    fGoodParticleTree->SetBranchAddress("BPROF4_Px",&fGoodBPROF4_Px);
    fGoodParticleTree->SetBranchAddress("BPROF4_Py",&fGoodBPROF4_Py);
    fGoodParticleTree->SetBranchAddress("BPROF4_Pz",&fGoodBPROF4_Pz);
    fGoodParticleTree->SetBranchAddress("BPROF4_PDGid",&fGoodBPROF4_PDGid);
    fGoodParticleTree->SetBranchAddress("BPROF4_EventID",&fGoodBPROF4_EventID);
    fGoodParticleTree->SetBranchAddress("BPROF4_TrackID",&fGoodBPROF4_TrackID);

    fGoodParticleTree->SetBranchAddress("BPROFEXT_x",&fGoodBPROFEXT_x);
    fGoodParticleTree->SetBranchAddress("BPROFEXT_y",&fGoodBPROFEXT_y);
    fGoodParticleTree->SetBranchAddress("BPROFEXT_z",&fGoodBPROFEXT_z);
    fGoodParticleTree->SetBranchAddress("BPROFEXT_t",&fGoodBPROFEXT_t);
    fGoodParticleTree->SetBranchAddress("BPROFEXT_Px",&fGoodBPROFEXT_Px);
    fGoodParticleTree->SetBranchAddress("BPROFEXT_Py",&fGoodBPROFEXT_Py);
    fGoodParticleTree->SetBranchAddress("BPROFEXT_Pz",&fGoodBPROFEXT_Pz);
    fGoodParticleTree->SetBranchAddress("BPROFEXT_PDGid",&fGoodBPROFEXT_PDGid);
    fGoodParticleTree->SetBranchAddress("BPROFEXT_EventID",&fGoodBPROFEXT_EventID);
    fGoodParticleTree->SetBranchAddress("BPROFEXT_TrackID",&fGoodBPROFEXT_TrackID);
         
    // Calculate the number of events to overlay
    CalculateNOverlays();
    
    // Now we need to fill the particle map
    FillParticleMaps();

    
    if( fSaveRecoTree ){
      fRecoTree = tfs->make<TTree>("tree", ""); 
      fRecoTree->Branch( "XBPF697_p", &fXBPF697_p );
      fRecoTree->Branch( "XBPF701_p", &fXBPF701_p );
      fRecoTree->Branch( "XBPF702_p", &fXBPF702_p );

      fRecoTree->Branch( "XBPF697_f", &fXBPF697_f );
      fRecoTree->Branch( "XBPF701_f", &fXBPF701_f );
      fRecoTree->Branch( "XBPF702_f", &fXBPF702_f );

      fRecoTree->Branch( "XBPF697_x", &fXBPF697_x );
      fRecoTree->Branch( "XBPF701_x", &fXBPF701_x );
      fRecoTree->Branch( "XBPF702_x", &fXBPF702_x );

      fRecoTree->Branch( "XBPF697_rx", &fXBPF697_rx );
      fRecoTree->Branch( "XBPF701_rx", &fXBPF701_rx );
      fRecoTree->Branch( "XBPF702_rx", &fXBPF702_rx );

      fRecoTree->Branch( "Reco_p",  &fReco_p );
      fRecoTree->Branch( "Reco_tof",  &fReco_tof );
      fRecoTree->Branch( "NP04front_p",  &fNP04front_p );

      fRecoTree->Branch( "NP04_PDG", &fNP04_PDG );

      fRecoTree->Branch( "XBPF707_x", &fXBPF707_x );
      fRecoTree->Branch( "XBPF707_f", &fXBPF707_f );
      fRecoTree->Branch( "XBPF708_y", &fXBPF708_y );
      fRecoTree->Branch( "XBPF708_f", &fXBPF708_f );
      fRecoTree->Branch( "XBPF716_x", &fXBPF716_x );
      fRecoTree->Branch( "XBPF716_f", &fXBPF716_f );
      fRecoTree->Branch( "XBPF717_y", &fXBPF717_y );
      fRecoTree->Branch( "XBPF717_f", &fXBPF717_f );

      fRecoTree->Branch( "XBPF707_rx", &fXBPF707_rx );
      fRecoTree->Branch( "XBPF708_ry", &fXBPF708_ry );
      fRecoTree->Branch( "XBPF716_rx", &fXBPF716_rx );
      fRecoTree->Branch( "XBPF717_ry", &fXBPF717_ry );

      fRecoTree->Branch( "TrueFront_x", &fTrueFront_x );
      fRecoTree->Branch( "TrueFront_y", &fTrueFront_y );
      fRecoTree->Branch( "TrueFront_z", &fTrueFront_z );

      fRecoTree->Branch( "TrueFront_Px", &fTrueFront_Px );
      fRecoTree->Branch( "TrueFront_Py", &fTrueFront_Py );
      fRecoTree->Branch( "TrueFront_Pz", &fTrueFront_Pz );

      fRecoTree->Branch( "RecoFront_x", &fRecoFront_x );
      fRecoTree->Branch( "RecoFront_y", &fRecoFront_y );
      fRecoTree->Branch( "RecoFront_z", &fRecoFront_z );
      fRecoTree->Branch("TrueID", &fTrueID);

    }
}

//----------------------------------------------------------------------------------------

void evgen::ProtoDUNEBeam::beginRun(art::Run& run)
{
    // Grab the geometry object to see what geometry we are using
    art::ServiceHandle<geo::Geometry> geo;
    std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetectorName()));
    run.put(std::move(runcol));
}

//--------------------------------------------------------------------------------------------

void evgen::ProtoDUNEBeam::endJob(){
    fInputFile->Close();
}


//--------------------------------------------------------------------------------------------
void evgen::ProtoDUNEBeam::produce(art::Event & e)
{
    
    // Define the truth collection for this event.
    auto truthcol = std::make_unique< std::vector<simb::MCTruth> >();
    
    //------------ Added by Caroline for beam simulation storage----------------------------
    std::unique_ptr<std::vector<sim::ProtoDUNEbeamsim>> beamsimcol (new std::vector<sim::ProtoDUNEbeamsim>);
    std::unique_ptr<std::vector<beam::ProtoDUNEBeamEvent> > beamData(new std::vector<beam::ProtoDUNEBeamEvent>);
    //std::unique_ptr<art::Assns<sim::ProtoDUNEbeamsim, simb::MCTruth> > beamsimassn (new art::Assns<sim::ProtoDUNEbeamsim, simb::MCTruth>);
    simb::MCTruth truth;
    beam::ProtoDUNEBeamEvent beamEvent;
    
    // Fill the MCTruth object
    GenerateTrueEvent(truth, (*beamsimcol), beamEvent ); // Add a reference for ProtoDUNEBeamEvent
    
    // Add the MCTruth to the vector
    truthcol->push_back(truth);
    
    //Make the assn                                                                                                  
    // util::CreateAssn(*this, e, *beamsimcol, truth, *beamsimassn);
    // Finally, add the MCTruth to the event
    e.put(std::move(truthcol));
    e.put(std::move(beamsimcol));

    beamData->push_back( beamEvent );
    e.put( std::move( beamData ) );

    //puts the vector object on to each event
    // We have made our event, increment the event number.
    ++fEventNumber;
}
//--------------------------------------------------------------------------------------

// Fill the particle maps using the input files. This links the events of interest
// to the entry number in fAllParticlesTree.
void evgen::ProtoDUNEBeam::FillParticleMaps(){
    
    // First off, loop over the good particles tree.
    int goodEventCounter = 0;
    for(int i = 0; i < fGoodParticleTree->GetEntries(); ++i){
        // If we want to skip some events, make sure we don't bother reading them in.
        if(fStartEvent > goodEventCounter){
            ++goodEventCounter;
            continue;
        }
        else{
            ++goodEventCounter;
        }
        
        fGoodParticleTree->GetEntry(i);

//        std::cout << "Tree entry " << i << " corresponds to event " << fGoodBPROFEXT_EventID << std::endl;

        // Make sure we didn't have two good particles in one event
        if(std::find(fGoodEventList.begin(),fGoodEventList.end(),(int)fGoodBPROFEXT_EventID)!=fGoodEventList.end()) continue;
 
        // NEW APPROACH - construct a ProtoFullSpill object
        ProtoFullSpill newSpill(fGoodNP04front_EventID,fGoodNP04front_TrackID,fGoodNP04front_t,i);
        fAllSpills.push_back(newSpill); 

        fGoodEventList.push_back(fGoodBPROFEXT_EventID);

    }
    
    // Print a message in case a user starts thinking something has broken.
    mf::LogInfo("ProtoDUNEBeam") << "About to loop over the beam simulation tree, this could take some time.";
    
    // Now we need to loop over the main particle tree
    for(int i = 0; i < fAllParticlesTree->GetEntries(); ++i){
        fAllParticlesTree->GetEntry(i);
        
        if (i%100000==0) std::cout << "Looking at entry " << i << std::endl;
        
        int event = int(fAllEventID);
       
        // Look at which good events this should be overlaid with
        std::vector<int> goodEventList = GetAllOverlays(event,fOverlays);

        unsigned int nMatches = 0;
        for(auto &spill : fAllSpills){
          // Stop looking if we have found all of our matches
          if(nMatches == goodEventList.size()) break;

          // Is this a spill of interest?
          if(std::find(goodEventList.begin(),goodEventList.end(),spill.fGoodEvent) == goodEventList.end()) continue;
          
          // Yes, so add this track to the spill
          spill.fAllSpillTracks[event].push_back(i);
          ++nMatches;
        }
        
    } // End loop over the main tree.
    
    mf::LogInfo("ProtoDUNEBeam") << "Found " << fGoodEventList.size() << " good events containing " << goodEventCounter << " good particles.";
    mf::LogInfo("ProtoDUNEBeam") << "Built " << fAllSpills.size() << " beam spills.";
    mf::LogInfo("ProtoDUNEBeam") << "All maps built, beginning event generation.";
    
}



//-------------------------------------------------------------------------------------------------
//modified by Caroline for beam sim storage
void evgen::ProtoDUNEBeam::GenerateTrueEvent(simb::MCTruth &mcTruth, std::vector<sim::ProtoDUNEbeamsim> &beamsimcol, beam::ProtoDUNEBeamEvent & beamEvent){
    //  std::unique_ptr<std::vector<sim::ProtoDUNEbeamsim>> beamsimcol
    // Check we haven't exceeded the length of the input tree
    if(fEventNumber >= (int)fGoodEventList.size()){
        throw cet::exception("ProtoDUNEBeam") << "Requested entry " << fEventNumber
        << " but tree only has entries 0 to "
        << fGoodEventList.size() - 1 << std::endl;
    } //end of if statement

    size_t nBeamEvents = 0;
    
    // Get the list of entries for the current event
//    int beamEvent = fGoodEventList[fCurrentGoodEvent];
    
    // A single particle seems the most accurate description.
    mcTruth.SetOrigin(simb::kSingleParticle);
   
    // NEW APPROACH
    ProtoFullSpill spill = fAllSpills[fCurrentGoodEvent];
//    std::cout << "This spill has " << spill.fAllSpillTracks.size() << " contributing events" << std::endl;

    // Find the entries that we are interested in.
        //std::cout << "Finding all particles associated with good particle event " << std::endl;
    for(auto const & event : spill.fAllSpillTracks){

        //std::cout << " - This event has " << event.second.size() << " contributing tracks" << std::endl;

        // Is this the event we would have triggered on?
        bool trigEvent = (event.first == spill.fGoodEvent);
        float baseTime;
        if(trigEvent){
            // Set the base time for the triggered event equal to the negative of the good particle time.
            // This will be corrected later on to set the time to zero, but keep time offsets within the event.
            baseTime = -1.0 * spill.fGoodTime;
        }
        else{
            // Get a random time from -fReadoutWindow to +fReadoutWindow in ns (fReadoutWindow value is in ms).
            baseTime = (fFlatRnd.fire() - 0.5)*2.0*(fReadoutWindow*1000.*1000.);
        }
        for(auto const t : event.second){
            // Get the entry from the tree for this event and track.
            fAllParticlesTree->GetEntry(t);
            
            //std::cout << fAllEventID << ", " << fAllTrackID << ", " << spill.fGoodEvent << ", " << spill.fGoodTrack << std::endl;

            // Convert the pdgCode to an int
            int intPDG = (int)fPDG;
            // We need to ignore nuclei for now...
            if(intPDG > 100000){ 
              //std::cout << "Skipping nuc" << std::endl;
              continue;
            }

            // Check to see if this should be a primary beam particle (good particle) or beam background
            std::string process="primaryBackground";

            TLorentzVector pos;
            TLorentzVector mom;
            // If this track is a "good particle", use the usual "primary" tag
            /*
            std::cout << "Trig: " << trigEvent << std::endl;
            std::cout << "SpillGoodEvent: " << spill.fGoodEvent << std::endl;
            std::cout << "SpillTrack: " << spill.fGoodTrack << std::endl;
            std::cout << "AllTrack: " << (int)fAllTrackID << std::endl;
            std::cout << "t: " << t << std::endl;
            */
            if(trigEvent && (spill.fGoodTrack == (int)fAllTrackID)){
              process="primary";
              // We also need to build the momentum vector using the correct good particle information
              fGoodParticleTree->GetEntry(spill.fGoodIndex);
              pos = ConvertCoordinates(fGoodNP04front_x/10.,fGoodNP04front_y/10.,fGoodNP04front_z/10.,baseTime + fGoodNP04front_t);
              mom = MakeMomentumVector(fGoodNP04front_Px/1000.,fGoodNP04front_Py/1000.,fGoodNP04front_Pz/1000.,(int)fGoodNP04front_PDGid,true);
              if (fSaveRecoTree) {
                fTrueFront_x = pos.X();
                fTrueFront_y = pos.Y();
                fTrueFront_z = pos.Z();

                fTrueFront_Px = mom.X();
                fTrueFront_Py = mom.Y();
                fTrueFront_Pz = mom.Z();
              }
              
              SetBeamEvent(beamEvent);
              ++nBeamEvents;
            }
            else{
              // We just need to shift our background particles upstream to BPROFEXT so they will hit the CRTs
              TVector3 tempPos = GetBackgroundPosition(fX,fY,fZ,fPx,fPy,fPz); 
              // At this step the position and momentum matches the GoodPartcle coordinates so apply the same functions
              pos = ConvertCoordinates(tempPos.X()/10.,tempPos.Y()/10.,tempPos.Z()/10.,baseTime+fEntryT);
              mom = MakeMomentumVector(fPx/1000.,fPy/1000.,fPz/1000.,intPDG,false);
              fGoodParticleTree->GetEntry(spill.fGoodIndex);
              std::cout << "SpillGoodTrack: " << spill.fGoodTrack << " " << (int)fGoodNP04front_PDGid << std::endl;
              std::cout << "ParentID: " << (int)fAllParentID << " " << (int)fPDG << std::endl;
              std::cout << "TrackID: " << (int)fAllTrackID<< std::endl;
              std::cout << "GoodEvent: " << spill.fGoodEvent << std::endl;
//              if(fabs(intPDG) == 13){
//                std::cout << "Found a " << process << " muon at time = " << baseTime+fEntryT << ":" << std::endl;
//                std::cout << fX/10. << ", " << fY/10. << ", " << fZ/10. << std::endl;
//                pos.Print();
//                mom.Vect().Unit().Print();
//              } 
            }

//            std::cout << "Information for particle " << intPDG << " with process " << process << std::endl;
//            pos.Print();
//            mom.Print();            
//            mom.Vect().Unit().Print();            
            
            // Track ID needs to be negative for primaries
            int trackID = -1*(mcTruth.NParticles() + 1); //g4trkid in larsoft
            if (fSaveRecoTree)
              fTrueID = trackID;
            
            // Create the particle and add the starting position and momentum
            //std::cout << "Adding particle with process " << process
            //          << " and PDG " << intPDG << std::endl;
            simb::MCParticle newParticle(trackID,intPDG,process);
            newParticle.AddTrajectoryPoint(pos,mom);
            
            // Add the MCParticle to the MCTruth for the event.
            mcTruth.Add(newParticle);

            // We want to save extra information from the beam monitors for the Good Particle
            if(trigEvent && (spill.fGoodTrack == (int)fAllTrackID)){

              fGoodParticleTree->GetEntry(spill.fGoodIndex);

              sim::ProtoDUNEBeamInstrument tof1("TOF1",fGoodTOF1_x,fGoodTOF1_y,fGoodTOF1_z,fGoodTOF1_t,fGoodTOF1_Px,fGoodTOF1_Py,fGoodTOF1_Pz,fGoodTOF1_PDGid,fGoodTOF1_EventID,fGoodTOF1_TrackID,fT_Resolution);  
              sim::ProtoDUNEBeamInstrument trig2("TRIG2",fGoodTRIG2_x,fGoodTRIG2_y,fGoodTRIG2_z,fGoodTRIG2_t,fGoodTRIG2_Px,fGoodTRIG2_Py,fGoodTRIG2_Pz,fGoodTRIG2_PDGid,fGoodTRIG2_EventID,fGoodTRIG2_TrackID,fT_Resolution);  

              // For BPROF4 we want to rotate the coordinates into the detector frame
              TLorentzVector bprof4Pos = ConvertBeamMonitorCoordinates(fGoodBPROF4_x,fGoodBPROF4_y,fGoodBPROF4_z,fGoodBPROF4_t,fBPROF4Pos);
              TVector3 bprof4Mom = ConvertBeamMonitorMomentumVec(fGoodBPROF4_Px,fGoodBPROF4_Py,fGoodBPROF4_Pz);
                                sim::ProtoDUNEBeamInstrument bprof4("BPROF4",bprof4Pos.X(),bprof4Pos.Y(),bprof4Pos.Z(),fGoodBPROF4_t,bprof4Mom.X(),bprof4Mom.Y(),bprof4Mom.Z(),fGoodBPROF4_PDGid,fGoodBPROF4_EventID,fGoodBPROF4_TrackID,fPos_Resolution);

              // Same for BPROFEXT
              TLorentzVector bprofextPos = ConvertBeamMonitorCoordinates(fGoodBPROFEXT_x,fGoodBPROFEXT_y,fGoodBPROFEXT_z,fGoodBPROFEXT_t,fBPROFEXTPos);
              TVector3 bprofextMom = ConvertBeamMonitorMomentumVec(fGoodBPROFEXT_Px,fGoodBPROFEXT_Py,fGoodBPROFEXT_Pz);
                                sim::ProtoDUNEBeamInstrument bprofext("BPROFEXT",bprofextPos.X(),bprofextPos.Y(),bprofextPos.Z(),fGoodBPROFEXT_t,bprofextMom.X(),bprofextMom.Y(),bprofextMom.Z(),fGoodBPROFEXT_PDGid,fGoodBPROFEXT_EventID,fGoodBPROFEXT_TrackID,fPos_Resolution);

//              std::cout << "Predicted detector position" << std::endl;
//              TVector3 predDir = (bprof4Pos.Vect()-bprofextPos.Vect()).Unit();
//              float projDist = (fNP04frontPos - fBPROF4Pos) + fabs(fBeamZ*10./fBMBasisZ.Z());
//              (bprof4Pos.Vect() + projDist*predDir).Print();

                        sim::ProtoDUNEBeamInstrument trig1("TRIG1",fGoodTRIG1_x,fGoodTRIG1_y,fGoodTRIG1_z,fGoodTRIG1_t,fGoodTRIG1_Px,fGoodTRIG1_Py,fGoodTRIG1_Pz,fGoodTRIG1_PDGid,fGoodTRIG1_EventID,fGoodTRIG1_TrackID,fT_Resolution);
                        sim::ProtoDUNEBeamInstrument bprof3("BPROF3",fGoodBPROF3_x,fGoodBPROF3_y,fGoodBPROF3_z,fGoodBPROF3_t,fGoodBPROF3_Px,fGoodBPROF3_Py,fGoodBPROF3_Pz,fGoodBPROF3_PDGid,fGoodBPROF3_EventID,fGoodBPROF3_TrackID,fPos_Resolution);
                        sim::ProtoDUNEBeamInstrument bprof2("BPROF2",fGoodBPROF2_x,fGoodBPROF2_y,fGoodBPROF2_z,fGoodBPROF2_t,fGoodBPROF2_Px,fGoodBPROF2_Py,fGoodBPROF2_Pz,fGoodBPROF2_PDGid,fGoodBPROF2_EventID,fGoodBPROF2_TrackID,fPos_Resolution);
                        sim::ProtoDUNEBeamInstrument bprof1("BPROF1",fGoodBPROF1_x,fGoodBPROF1_y,fGoodBPROF1_z,fGoodBPROF1_t,fGoodBPROF1_Px,fGoodBPROF1_Py,fGoodBPROF1_Pz,fGoodBPROF1_PDGid,fGoodBPROF1_EventID,fGoodBPROF1_TrackID,fPos_Resolution);

// Adding Cherenkovs with same variables as BPROFEXT except for their response
sim::ProtoDUNEBeamInstrument cherenkov1("CHERENKOV1",fGoodBPROFEXT_x,fGoodBPROFEXT_y,fGoodBPROFEXT_z,fGoodBPROFEXT_t,fGoodBPROFEXT_Px,fGoodBPROFEXT_Py,fGoodBPROFEXT_Pz,fGoodBPROFEXT_PDGid,fGoodBPROFEXT_EventID,fGoodBPROFEXT_TrackID,fCh_Efficiency);
sim::ProtoDUNEBeamInstrument cherenkov2("CHERENKOV2",fGoodBPROFEXT_x,fGoodBPROFEXT_y,fGoodBPROFEXT_z,fGoodBPROFEXT_t,fGoodBPROFEXT_Px,fGoodBPROFEXT_Py,fGoodBPROFEXT_Pz,fGoodBPROFEXT_PDGid,fGoodBPROFEXT_EventID,fGoodBPROFEXT_TrackID,fCh_Efficiency);


                        sim::ProtoDUNEbeamsim temp; 
                        temp.AddInstrument(tof1);
                        temp.AddInstrument(trig2);
                        temp.AddInstrument(bprof4);
                        temp.AddInstrument(bprofext);
                        temp.AddInstrument(trig1);
                        temp.AddInstrument(bprof3);
                        temp.AddInstrument(bprof2);
                        temp.AddInstrument(bprof1);
                        temp.AddInstrument(cherenkov1);
                        temp.AddInstrument(cherenkov2);

//                        std::cout << "ProtoDUNEbeamsim object has " << temp.NInstruments() << " beam instruments" << std::endl;
//std::cout << "TOF1 resolution: " << fT_Resolution << std::endl;

                        beamsimcol.push_back(temp);
//                        std::cout<< beamsimcol.size() << std::endl;
                        // std::cout<<" test value beam profile monitor: TTREE   "<<fGoodBPROF4_x<<std::endl;
//                        std::cout<<"From TTree TRIG2_TRACKID: "<<fGoodTRIG2_TrackID<<std::endl;
//                        std::cout<<"From ProtoDUNEBeamInstrument: "<<tof1.GetT()<<std::endl;
//                        std::cout<<"From ProtoDUNEBeamInstrument: "<<tof1.GetSmearedVar1()<<std::endl;


                        // std::cout<<"the testing for beam profile monitor information:  "<<fGoodBPROF4_z<<std::endl;
                        
                        
//                        std::cout<< "From the data product:  TRIG2TRACKID:   "<<temp.get_TRIG2_TrackID()<<std::endl;
                        //check the last index of the vector
                        sim::ProtoDUNEbeamsim lastelement = beamsimcol.back();
                        
//                        std::cout<<"From the vector TRIG2_TRACKID: "<<lastelement.get_TRIG2_TrackID()<<std::endl;
                        
                        
                        //Make the assn                                                                                                                                  
                        //util::CreateAssn(*this, e, *beamsimcol, newParticle, *beamsimassn)
                        

            } // End beam instrumentation section
            
        } // End loop over interesting tracks for each event
    } // End loop over the vector of interesting events

    mf::LogInfo("ProtoDUNEBeam") << "Got " << nBeamEvents << " beam events";
    mf::LogInfo("ProtoDUNEBeam") << "Created event with " << mcTruth.NParticles() << " particles.";
    
    // Move on the good event iterator
    ++fCurrentGoodEvent;
}

//---------------------------------------------------------------------------------------

// Function written in similar way as "openDBs()" in CORSIKAGen_module.cc
void evgen::ProtoDUNEBeam::OpenInputFile()
{
    // Setup ifdh object
    if (!fIFDH)
    {
        fIFDH = new ifdh_ns::ifdh;
    }
    
    const char* ifdh_debug_env = std::getenv("IFDH_DEBUG_LEVEL");
    if ( ifdh_debug_env )
    {
        mf::LogInfo("ProtoDUNEBeam") << "IFDH_DEBUG_LEVEL: " << ifdh_debug_env<<"\n";
        fIFDH->set_debug(ifdh_debug_env);
    }
    
    std::string path(gSystem->DirName(fFileName.c_str()));
    std::string pattern(gSystem->BaseName(fFileName.c_str()));
    
    auto flist = fIFDH->findMatchingFiles(path,pattern);
    if (flist.empty())
    {
        struct stat buffer;
        if (stat(fFileName.c_str(), &buffer) != 0)
        {
            throw cet::exception("ProtoDUNEBeam") << "No files returned for path:pattern: "<<path<<":"<<pattern<<std::endl;
        }
        else
        {
            mf::LogInfo("ProtoDUNEBeam") << "For "<< fFileName <<"\n";
        }
    }
    else
    {
        std::pair<std::string, long> f = flist.front();
        
        mf::LogInfo("ProtoDUNEBeam") << "For "<< fFileName <<"\n";
        
        // Do the fetching, store local filepaths in locallist
        
        mf::LogInfo("ProtoDUNEBeam")
        << "Fetching: " << f.first << " " << f.second <<"\n";
        std::string fetchedfile(fIFDH->fetchInput(f.first));
        MF_LOG_DEBUG("ProtoDUNEBeam") << " Fetched; local path: " << fetchedfile;
        
        fFileName = fetchedfile;
    }
}


//----------------------------------------------------------------------------------

TLorentzVector evgen::ProtoDUNEBeam::ConvertCoordinates(float x, float y, float z, float t){
    
    float finalX = x + fBeamX;
    float finalY = y + fBeamY;
//    float finalZ = (z - z) + fBeamZ; // Just use the z position
    float finalZ = z + fBeamZ; // Just use the z position
    
    TLorentzVector newPos(finalX,finalY,finalZ,t);
    return newPos;
}

//--------------------------------------------------------------------------------

TLorentzVector evgen::ProtoDUNEBeam::MakeMomentumVector(float px, float py, float pz, int pdg, bool shifts){
    
    //float rotationXZ = fRotateXZ;
    //float rotationYZ = fRotateYZ;
    
    // Make the momentum vector and rotate it
    TVector3 momVec(px,py,pz);

    if(shifts){
      // Shift theta and phi first, and only rotate if we should
      if(fabs(fBeamThetaShift) > 1.e-10 && fabs(fBeamPhiShift) > 1.e-10){
        //std::cout << "Shifted momentum vector from " << momVec.Theta() << ", " << momVec.Phi() << " ";
        momVec.SetTheta(momVec.Theta() + fBeamThetaShift);
        momVec.SetPhi(momVec.Phi() + fBeamPhiShift);
        //std::cout << "to " << momVec.Theta() << ", " << momVec.Phi() << std::endl;
      }
      else{
        momVec.RotateY(fRotateXZ/*rotationXZ*/ * TMath::Pi() / 180.);
        momVec.RotateX(fRotateYZ/*rotationYZ*/ * TMath::Pi() / 180.);
      }
    }

    // Find the particle mass so we can form the energy
    const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
    const TParticlePDG* definition = databasePDG->GetParticle(pdg);
    float mass = definition->Mass();
    
    float energy = sqrt(mass*mass + momVec.Mag2());
    
    TLorentzVector newMom(momVec,energy);
    return newMom;
}

TLorentzVector evgen::ProtoDUNEBeam::MakeMomentumVector(const TVector3 &mom, int pdg, bool shifts){

  return MakeMomentumVector(mom.X(),mom.Y(),mom.Z(),pdg,shifts);

}

//-----------------------------------------------------------------------------

void evgen::ProtoDUNEBeam::CalculateNOverlays(){
    
    // The number of events to overlay is as follows:
    // N = Intensity * 2.0 * ReadoutWindow / BeamSpillLength
    fOverlays = fIntensity * (2.0 * fReadoutWindow / 1000.) / fBeamSpillLength;
    std::cout << "Number of overlays = " << fOverlays << std::endl;   
}


//-------------------------------------------------------------------------------
int evgen::ProtoDUNEBeam::IsOverlayEvent(int event, int nOverlay){
    
    // Check if this event lies within nOverlay/2 of each
    for(auto const e : fGoodEventList){
        if(fabs(event - e) < nOverlay/2){
            return e;
        }
    }
    return -1;
}


//---------------------------------------------------------------------------------
std::vector<int> evgen::ProtoDUNEBeam::GetAllOverlays(int event, int nOverlay){
    
    std::vector<int> nMatches;
    for(auto const e : fGoodEventList){
        if(fabs(event - e) < nOverlay/2){
            nMatches.push_back(e);
        }
    }
    return nMatches;
    
}
//----------------------------------------------------------------------------------

// We need to rotate the beam monitor coordinates into the detector frame (matching NP04front)
// This means they can later be treated in the same way as the standard NP04front positions
TLorentzVector evgen::ProtoDUNEBeam::ConvertBeamMonitorCoordinates(float x, float y, float z, float t, float zOffset){

  float off = fNP04frontPos - zOffset;

  TLorentzVector old(x,y,z,t);

  // Convert the coordinates using the rotated basis vectors
  float newX = x*fBMBasisX.X() + y*fBMBasisY.X() + (z-zOffset)*fBMBasisZ.X() + off*fabs(fBMBasisZ.X());
  float newY = x*fBMBasisX.Y() + y*fBMBasisY.Y() + (z-zOffset)*fBMBasisZ.Y() + off*fabs(fBMBasisZ.Y());
  float newZ = x*fBMBasisX.Z() + y*fBMBasisY.Z() + (z-zOffset) - off*fabs(fBMBasisZ.Z());

  // Account for the small differences between NP04front and the detector coordinates
  newX += fBeamX*10.;
  newY += fBeamY*10.;
  newZ += fBeamZ*10.;

  // Make our new beam monitor position in the detector coordinate system
  TLorentzVector result(newX,newY,newZ,t);

//  std::cout << "Coordinate transform..." << std::endl;
//  old.Print();
//  result.Print();

  return result;
}

TVector3 evgen::ProtoDUNEBeam::ConvertProfCoordinates(double x, double y, double z, double zOffset){
  double off = fNP04frontPos - zOffset;

  TVector3 old(x,y,z);

  double newX = x*fBMBasisX.X() + y*fBMBasisY.X() + /*(z-zOffset)*fBMBasisZ.X()*/ + off*fabs(fBMBasisZ.X());
  double newY = x*fBMBasisX.Y() + y*fBMBasisY.Y() + /*(z-zOffset)*fBMBasisZ.Y()*/ + off*fabs(fBMBasisZ.Y());
  double newZ = x*fBMBasisX.Z() + y*fBMBasisY.Z() + /*(z-zOffset)              */ - off*fabs(fBMBasisZ.Z());

  newX += fBeamX*10.;
  newY += fBeamY*10.;
  newZ += fBeamZ*10.;

  TVector3 result(newX/10., newY/10., newZ/10.);
  return result;
}

TVector3 evgen::ProtoDUNEBeam::ConvertBeamMonitorMomentumVec(float px, float py, float pz){

  TVector3 newMom(px,py,pz);
//  std::cout << "Momentum transform..." << std::endl;
//  newMom.Unit().Print();
  RotateMonitorVector(newMom);
//  newMom.Unit().Print();
  return newMom;
}

void evgen::ProtoDUNEBeam::BeamMonitorBasisVectors(){

  fBMBasisX = TVector3(1.,0.,0.);
  fBMBasisY = TVector3(0.,1.,0.);
  fBMBasisZ = TVector3(0.,0.,1.);
  RotateMonitorVector(fBMBasisX);
  RotateMonitorVector(fBMBasisY);
  RotateMonitorVector(fBMBasisZ);

}

void evgen::ProtoDUNEBeam::RotateMonitorVector(TVector3 &vec){

  // Note: reordering how these are done in order to keep the basis
  //       vectors of the monitors parallel to the ground. 
  vec.RotateX( fRotateMonitorYZ * TMath::Pi() / 180. );
  vec.RotateY( fRotateMonitorXZ * TMath::Pi() / 180. );

}

TVector3 evgen::ProtoDUNEBeam::GetBackgroundPosition(float x, float y, float z, float px, float py, float pz){

  TVector3 pos(x,y,z);
  TVector3 dir = TVector3(px,py,pz).Unit();

  // Want to move the position upstream by a distance equal to fNP04frontPos - fBPROFEXTPos
  // This length is in the beam direction frame unless we account for it
  float shiftLength = (fNP04frontPos - fBPROFEXTPos)/fBMBasisZ.Z();

  return pos - shiftLength*dir;
}


void evgen::ProtoDUNEBeam::SetBeamEvent(beam::ProtoDUNEBeamEvent & beamevt){
  //This will just use the class members
  
  beamevt.SetTOFs( std::vector<double>{ fGoodTRIG2_t - fGoodTOF1_t } );
  beamevt.SetTOFChans( std::vector<int>{ 0 } );
  beamevt.SetUpstreamTriggers( std::vector<size_t>{0} );
  beamevt.SetDownstreamTriggers( std::vector<size_t>{0} );
  beamevt.SetCalibrations( 0., 0., 0., 0. );
  beamevt.DecodeTOF();

  beamevt.SetMagnetCurrent( 0. );
  beamevt.SetTimingTrigger( 12 );

  beam::CKov dummy;
  dummy.trigger = 0;
  dummy.pressure = 0.;
  dummy.timeStamp = 0.;
  beamevt.SetCKov0( dummy );
  beamevt.SetCKov1( dummy );

  beamevt.SetActiveTrigger(0);
  beamevt.SetT0( std::make_pair(0.,0.) );

  beamevt.SetFBMTrigger( "XBPF022697", MakeFiberMonitor( fGoodBPROF1_x ) );
  beamevt.SetFBMTrigger( "XBPF022698", MakeFiberMonitor( fGoodBPROF1_y ) );
  beamevt.SetFBMTrigger( "XBPF022701", MakeFiberMonitor( fGoodBPROF2_x ) );
  beamevt.SetFBMTrigger( "XBPF022702", MakeFiberMonitor( fGoodBPROF3_x ) );

  beamevt.SetFBMTrigger( "XBPF022707", MakeFiberMonitor( fGoodBPROFEXT_x ) );
  beamevt.SetFBMTrigger( "XBPF022708", MakeFiberMonitor( fGoodBPROFEXT_y ) );

  beamevt.SetFBMTrigger( "XBPF022716", MakeFiberMonitor( fGoodBPROF4_x ) );
  beamevt.SetFBMTrigger( "XBPF022717", MakeFiberMonitor( fGoodBPROF4_y ) );

  MakeTracks( beamevt );
  MomentumSpectrometer( beamevt );


  if( fSaveRecoTree ){ 
    fReco_p = beamevt.GetRecoBeamMomentum(0);
    fReco_tof = beamevt.GetTOF();
    std::cout << "TOF: " << beamevt.GetTOFs()[0] << " " << beamevt.GetTOF() << std::endl;
    fNP04_PDG = fGoodNP04front_PDGid;

    fNP04front_p = sqrt( fGoodNP04front_Px * fGoodNP04front_Px 
                       + fGoodNP04front_Py * fGoodNP04front_Py 
                       + fGoodNP04front_Pz * fGoodNP04front_Pz );

    fXBPF697_p = sqrt( fGoodBPROF1_Px * fGoodBPROF1_Px 
                     + fGoodBPROF1_Py * fGoodBPROF1_Py 
                     + fGoodBPROF1_Pz * fGoodBPROF1_Pz );

    fXBPF701_p = sqrt( fGoodBPROF2_Px*fGoodBPROF2_Px 
                     + fGoodBPROF2_Py*fGoodBPROF2_Py 
                     + fGoodBPROF2_Pz*fGoodBPROF2_Pz );

    fXBPF702_p = sqrt( fGoodBPROF3_Px*fGoodBPROF3_Px 
                     + fGoodBPROF3_Py*fGoodBPROF3_Py 
                     + fGoodBPROF3_Pz*fGoodBPROF3_Pz );
    fXBPF697_x = fGoodBPROF1_x;
    fXBPF701_x = fGoodBPROF2_x;
    fXBPF702_x = fGoodBPROF3_x;

    fXBPF716_x = fGoodBPROF4_x;
    fXBPF717_y = fGoodBPROF4_y;
    fXBPF707_x = fGoodBPROFEXT_x;
    fXBPF708_y = fGoodBPROFEXT_y;

    
    fXBPF697_f = beamevt.GetFBM( "XBPF022697" ).active[0];
    fXBPF701_f = beamevt.GetFBM( "XBPF022701" ).active[0];
    fXBPF702_f = beamevt.GetFBM( "XBPF022702" ).active[0];

    fXBPF707_f = beamevt.GetFBM( "XBPF022707" ).active[0];
    fXBPF708_f = beamevt.GetFBM( "XBPF022708" ).active[0];
    fXBPF716_f = beamevt.GetFBM( "XBPF022716" ).active[0];
    fXBPF717_f = beamevt.GetFBM( "XBPF022717" ).active[0];

    fXBPF697_rx = GetPosition( fXBPF697_f ); 
    fXBPF701_rx = GetPosition( fXBPF701_f ); 
    fXBPF702_rx = GetPosition( fXBPF702_f ); 
                                         
    fXBPF716_rx = GetPosition( fXBPF716_f ); 
    fXBPF717_ry = GetPosition( fXBPF717_f ); 
    fXBPF707_rx = GetPosition( fXBPF707_f ); 
    fXBPF708_ry = GetPosition( fXBPF708_f ); 

    //fTrueFront_x = fGoodNP04front_x + fBeamX;
    //fTrueFront_y = fGoodNP04front_y + fBeamY;
    //fTrueFront_z = fGoodNP04front_z + fBeamZ;

    fRecoFront_x = beamevt.GetBeamTrack(0).End().X();
    fRecoFront_y = beamevt.GetBeamTrack(0).End().Y();
    fRecoFront_z = beamevt.GetBeamTrack(0).End().Z();

    fRecoTree->Fill();
  }
}

beam::FBM evgen::ProtoDUNEBeam::MakeFiberMonitor( float pos ){
  beam::FBM theFBM;

  //I should probably just make this into
  //a constructor for the FBM...
  theFBM.ID = -1;
  theFBM.glitch_mask = {};
  std::uninitialized_fill( std::begin(theFBM.fiberData), std::end(theFBM.fiberData), 0. );
  std::uninitialized_fill( std::begin(theFBM.timeData), std::end(theFBM.timeData), 0. );
  theFBM.timeStamp = 0.;

  short f = 96 -  short( floor(pos) ) - 1;
  theFBM.fibers[f] = 1;
  theFBM.active.push_back(f);
  theFBM.decoded = true;

  return theFBM; 
}

double evgen::ProtoDUNEBeam::GetPosition( short fiber ){
  return ((96 - fiber) - .5);
}

void evgen::ProtoDUNEBeam::MakeTracks( beam::ProtoDUNEBeamEvent & beamEvent ){
  
  //We should only have one active fiber at a time
  //
  //Might need to ask Leigh, etc. if it's possible
  //to have multiple particles going through at the
  //same time. In which case -- try to implement it
  
  short fx1 = beamEvent.GetFBM( "XBPF022707" ).active[0];
  short fy1 = beamEvent.GetFBM( "XBPF022708" ).active[0];

  double x1 = GetPosition( fx1 );
  double y1 = GetPosition( fy1 );

  TVector3 pos1 = ConvertProfCoordinates( x1, y1, 0., fBPROFEXTPos );

  short fx2 = beamEvent.GetFBM( "XBPF022716" ).active[0];
  short fy2 = beamEvent.GetFBM( "XBPF022717" ).active[0];

  double x2 = GetPosition( fx2 );
  double y2 = GetPosition( fy2 );

  TVector3 pos2 = ConvertProfCoordinates( x2, y2, 0., fBPROF4Pos );
 
  std::vector< TVector3 > thePoints = { pos1, pos2, ProjectToTPC( pos1, pos2 ) };
  std::vector< TVector3 > theMomenta = {
    ( pos2 - pos1 ).Unit(),
    ( pos2 - pos1 ).Unit(),
    ( pos2 - pos1 ).Unit()
  };

  beamEvent.AddBeamTrack(
    recob::Track(
      recob::TrackTrajectory(recob::tracking::convertCollToPoint( thePoints ),
                             recob::tracking::convertCollToVector( theMomenta ),
                             recob::Track::Flags_t( thePoints.size() ),
                             false ),
      0, -1., 0, recob::tracking::SMatrixSym55(), recob::tracking::SMatrixSym55(), 1 
    )
  );
    
}

TVector3 evgen::ProtoDUNEBeam::ProjectToTPC(TVector3 firstPoint, TVector3 secondPoint){
  TVector3 dR = (secondPoint - firstPoint);
  
  double deltaZ = -1.*secondPoint.Z();
  double deltaX = deltaZ * (dR.X() / dR.Z());
  double deltaY = deltaZ * (dR.Y() / dR.Z());

  TVector3 lastPoint = secondPoint + TVector3(deltaX, deltaY, deltaZ);
  return lastPoint;
}

void evgen::ProtoDUNEBeam::MomentumSpectrometer( beam::ProtoDUNEBeamEvent & beamEvent ){

  short f1 = beamEvent.GetFBM( "XBPF022697" ).active[0];
  short f2 = beamEvent.GetFBM( "XBPF022701" ).active[0];
  short f3 = beamEvent.GetFBM( "XBPF022702" ).active[0];

  double x1 = -1.e-3 * GetPosition( f1 );
  double x2 = -1.e-3 * GetPosition( f2 );
  double x3 = -1.e-3 * GetPosition( f3 );

  double cos_theta = MomentumCosTheta( x1, x2, x3 );
  double momentum = 299792458*fLB/(1.E9 * acos(cos_theta));
  beamEvent.AddRecoBeamMomentum( momentum );

}

double evgen::ProtoDUNEBeam::MomentumCosTheta(double x1, double x2, double x3){
  double a =  (x2*fL3 - x3*fL2)*cos(fBeamBend)/(fL3-fL2);

 
  double numTerm = (a - x1)*( (fL3 - fL2)*tan(fBeamBend) + (x3 - x2)*cos(fBeamBend) ) + fL1*( fL3 - fL2 );

  double denomTerm1, denomTerm2, denom;
  denomTerm1 = sqrt( fL1*fL1 + (a - x1)*(a - x1) );
  denomTerm2 = sqrt( TMath::Power( ( (fL3 - fL2)*tan(fBeamBend) + (x3 - x2)*cos(fBeamBend) ),2)
                   + TMath::Power( ( (fL3 - fL2) ),2) );
  denom = denomTerm1 * denomTerm2;

  double cosTheta = numTerm/denom;  
  return cosTheta;
}

DEFINE_ART_MODULE(evgen::ProtoDUNEBeam)