CosmicsGen_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   Generator for cosmic-rays
/// \author  brebel@fnal.gov
/// \date
////////////////////////////////////////////////////////////////////////

extern "C" {
#include <getopt.h>
}
#include <cmath>  // for fabs() etc
#include <memory>

// ROOT includes
#include "TH1.h"
#include "TH2.h"

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/SubRun.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvA includes

#include "nutools/EventGeneratorBase/evgenbase.h"
#include "nutools/EventGeneratorBase/CRY/CRYHelper.h"
#include "Geometry/Geometry.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "SummaryData/CosmicExposure.h"
#include "SummaryData/RunData.h"
#include "SummaryData/SubRunData.h"


namespace evgen {

  /// A module to produce cosmic ray events
  class CosmicsGen : public art::EDProducer {
  public:
    explicit CosmicsGen(fhicl::ParameterSet const& pset);
    virtual ~CosmicsGen();                        

    void produce(art::Event& evt);
    virtual void beginRun(art::Run& run);
    virtual void beginSubRun(art::SubRun& run);
    virtual void endSubRun(art::SubRun& run);

  private:

    /// Method to cut the events not going through the DetectorBigBox,
    /// which is defined as the box surrounding the DetectorEncosureBox
    /// by adding an additional length in each of the dimensions.
    /// The additional length is defined in detectrorbigbox.xml
    void DetectorBigBoxCut(simb::MCTruth& truth);
    
    
    /// Method to project the cosmic particles back to the surface for the near detector.
    /// CRYHelper is setup such that it projects CRY particles, which generated
    /// at the surface up to the end of the World volume.
    /// This is because ndos, for instance, is above the surface.
    /// CRYHelper is standard (in the EventGeneratorBase), used by other experiments.
    /// So, we don't want to modify it.
    /// But we do project cosmics back to the surface for the near detector
    /// (maybe should do the same for the far???)
    void ProjectCosmicsToSurface(simb::MCTruth& truth);
    
    /// Method to project muons generated at the near detector to the DetectorBigBox.
    /// This is to get rid of the muon propagation through the rock, since
    /// the secondary particles would probably not reach the detector.
    /// Also we account for the muon's energy loss.
    /// After the muon reaches the DetectorBigBox, we start propagating it through Geant.
    void ProjectMuonsToDetectorBigBox(simb::MCTruth& truth);

    /// Method to take the MCTruth object and beak it into pieces for each individual Cosmic Ray
    void SeparateMCTruth(simb::MCTruth& truth, std::vector<simb::MCTruth>* truthcol);

    /// Method to produce and enhanced sample with a richer population of high energy particles
    /// with PDG set in fEnhancePDG. This is done at this step (before Geant) to minimize running 
    /// time. It only decides to keep or throw out complete events based on their most interesting
    /// particle so as to not interfere with the CRYHelper process.
    /// Currently set for neutrons, photons or protons independently.
    void Enhance(simb::MCTruth& truth, bool* skipThisSet, double random_number);
    
    bool isIntersectTheBox(const double xyz[],
                           const double dxyz[],
                           double xlo, double xhi,
                           double ylo, double yhi,
                           double zlo, double zhi);
    
    void ProjectToSurface (const double xyz[],
                           const double dxyz[],
                           double xlo, double xhi,
                           double ylo, double yhi,
                           double zlo, double zhi,
                           const double surface_y,
                           double xyzout[]);
    
    void ProjectToBox      (const double xyz[],
                            const double dxyz[],
                            double xlo, double xhi,
                            double ylo, double yhi,
                            double zlo, double zhi,
                            double xyzout[], double& dlout);
    
  private:

    CLHEP::HepRandomEngine & fEngine;
    evgb::CRYHelper* fCRYHelp;            ///< CRY generator object                     
    bool             fIsBigBoxUsed;     ///< Do we need to use the BigBox cut?
                                        ///< The cosmic rays must go through the DetectorBigBox.
                                        ///< Otherwise, don't store them
    bool             fBreakMCTruth;
    double           fExposure;         ///< Livetime this subrun, in seconds
    bool             fGenSingleEvents;  ///< generate one cosmic per record
    double           fProjectedZOffset; ///< How far above the BigBox ND muons should be started
    bool             fEnhance;          ///< Flag to enhace high energy particles
    int              fEnhancePDG;       ///< PDG of enhanced particles (neutron, proton or photon)
    int              fCycle;            ///< MC cycle generation number.
  };
}

////////////////////////////////////////////////////////////////////////

namespace evgen{

  //____________________________________________________________________________
  CosmicsGen::CosmicsGen(fhicl::ParameterSet const& pset)
    : EDProducer(pset)
    , fEngine(createEngine(pset.get< int >("Seed", evgb::GetRandomNumberSeed())))
    , fIsBigBoxUsed       (pset.get< bool   >("IsBigBoxUsed"))
    , fBreakMCTruth       (pset.get< bool   >("BreakMCTruth"))
    , fExposure           (0)
    , fGenSingleEvents    (pset.get< bool   >("GenSingleEvents"))
    , fProjectedZOffset   (pset.get< double >("ProjectedZOffset"))
    , fEnhance            (pset.get< bool   >("Enhance"))
    , fEnhancePDG         (pset.get< int    >("EnhancePDG"))
    , fCycle              (pset.get< int    >("Cycle", 0))
  {
    
    // Create a new CRYHelper
    fCRYHelp = new evgb::CRYHelper(pset, fEngine);
    
    // Producers
    produces< std::vector<simb::MCTruth> >();
    produces< sumdata::SubRunData, art::InSubRun >();
    produces< sumdata::RunData,    art::InRun    >();
    
    // This information is misleading when we're generating in normal mode, with
    // many cosmics per spill. CosmicExposureInfo gets it right. But for single
    // event mode we need it to do the accounting for overlays.
    // We also need accounting of exposure for hadron enhanced samples.
    if(fGenSingleEvents || fEnhance){
      produces< sumdata::CosmicExposure, art::InSubRun >();
    } 
 }

  //____________________________________________________________________________
  CosmicsGen::~CosmicsGen()
  {
    if(fCRYHelp) delete fCRYHelp;
  }

  //____________________________________________________________________________
  void CosmicsGen::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->DetId(),
                                                                geo->FileBaseName(),
                                                                geo->ExtractGDML()));
    
    run.put(std::move(runcol));

    return;
  }

  //____________________________________________________________________________
  void CosmicsGen::beginSubRun(art::SubRun& /*subrun*/)
  {
    fExposure = 0;
  }

  //____________________________________________________________________________
  void CosmicsGen::endSubRun(art::SubRun& subrun)
  {
    if(fGenSingleEvents || fEnhance){
      std::unique_ptr<sumdata::CosmicExposure> ce(new sumdata::CosmicExposure);
      ce->totlivetime = fExposure;
      subrun.put(std::move(ce));
    }
    
    // store the cycle information
    art::ServiceHandle<geo::Geometry> geo;
    std::unique_ptr< sumdata::SubRunData > sd(new sumdata::SubRunData(fCycle));
    subrun.put(std::move(sd));

  }

  //____________________________________________________________________________
  void CosmicsGen::produce(art::Event& evt)
  {
    std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);

    art::ServiceHandle<geo::Geometry> geo;
    double rantime = fEngine.flat();
 
    simb::MCTruth  truth;
    truth.SetOrigin(simb::kCosmicRay);

    double sampletime;
    double total_sampletime = 0.;

    // While skippedSet was implemented to regenerate the Truth object multiple
    // times for enhanced samples but is skipped for normal CRY.
    bool skippedSet = true;
    int total_skipped = -1;
    while ( skippedSet ){
      total_skipped++;
        
      rantime = fEngine.flat();
      sampletime = fCRYHelp->Sample(truth, geo->SurfaceY(), geo->DetLength(), 0, rantime);
        
      // 1. Project the cosmics rays down to the surface (ND only)
      this->ProjectCosmicsToSurface(truth);
      
      // 2. DetectorBigBox cut
      if(fIsBigBoxUsed) this->DetectorBigBoxCut(truth);
      
      // 3. ND only: Rewind cosmics back to the surface after adjustment by 2.
      this->ProjectCosmicsToSurface(truth);
      
      total_sampletime += sampletime;

      // * Option to produce an enhanced sample.
      if(!fEnhance) break;
      double rand_num = fEngine.flat();
      this->Enhance( truth, &skippedSet, rand_num);
      
    }

    // 4. ND only: Now project muons to the Big Box
    this->ProjectMuonsToDetectorBigBox(truth);
    
    //5. Break MCTruth into individual truth objects for each cosmic
    if(fBreakMCTruth){
      this->SeparateMCTruth(truth, truthcol.get());
    }
    else {
      truthcol->push_back(truth);
    }

    evt.put(std::move(truthcol));

    fExposure += total_sampletime;

    return;

  }

  ////////////////////////////////////////////////////////////////////////
  // Method to enhance high E Neutrons, Protons or Photons by deciding 
  // whether to keep or skip the current Truth based on a probability. 
  void CosmicsGen::Enhance(simb::MCTruth& truth, bool* skipThisSet, double random_number){

    // Parameters from the energy distribution fits docdb 13199
    // TO DO: Make these fcl parameters

    double nNeutron  = -5.18;
    double E0Neutron = 4;

    double nPhoton  = -2.63;
    double E0Photon = 3;

    double nProton  = -3.91;
    double E0Proton = 4;

    // Enhance with weight w = k(E/E0)^n 
    double n, E0;

    if( fEnhancePDG == 2112 ){
      n  = nNeutron;
      E0 = E0Neutron;
    }
    if( fEnhancePDG == 2212 ){
      n  = nProton;
      E0 = E0Proton;
    }
    if( fEnhancePDG == 22 ){
      n  = nPhoton;
      E0 = E0Photon;
    }

    double event_score = 0.001; // probability to keep evt, better for high E fEnhancePDGs
    double event_w = 1.;        // assign weight to relate these events to normal CRY stats

    int nparticles = truth.NParticles();
    for ( int part_Idx = 0; part_Idx < nparticles; ++part_Idx){
      simb::MCParticle particle = truth.GetParticle(part_Idx);
      int pdg = abs(particle.PdgCode());
      // Calculate event score by looking at interesting particles

      if ( pdg == fEnhancePDG ){
        const TLorentzVector& fourMomentum = particle.Momentum();
        double energy = fourMomentum.E();
        // Keep the largest particle score for the event
        double score = pow((E0/energy),n);
        if ( score > event_score ) event_score = score;
        // Keep all events with score >= 1
        if ( event_score > 1 ) break;
      }

    }// All particles in truth
    if ( event_score>1 ) event_score = 1;
    event_w = 1.0/event_score;

    simb::MCTruth  Enhancedtruth;
    Enhancedtruth.SetOrigin(simb::kCosmicRay);
    if ( random_number<event_score ) {
      *skipThisSet = false;
      int npart = truth.NParticles();
      for( int ipart = 0; ipart < npart; ++ipart ){
        simb::MCParticle particle = truth.GetParticle(ipart);

        simb::MCParticle part(particle.TrackId(),
                              particle.PdgCode(),
                              particle.Process(),
                              particle.Mother(),
                              particle.Mass(),
                              particle.StatusCode());
        part.AddTrajectoryPoint(particle.Position(),
                                particle.Momentum());
        part.SetWeight(event_w);

        Enhancedtruth.Add(part);
      }
    }
    // Enhancedtruth remains empty unless we keep the event.
    truth=Enhancedtruth;
  }

  ////////////////////////////////////////////////////////////////////////
  // Method to cut the events not going through the DetectorBigBox,
  // which is defined as the box surrounding the DetectorEncosureBox
  // by adding an additional length in each of the dimensions.
  // The additional length is defined in detectrorbigbox.xml
  void CosmicsGen::DetectorBigBoxCut(simb::MCTruth& truth){
    simb::MCTruth mctruth_new;
    
    double x1=0;
    double x2=0;
    double y1=0;
    double y2=0;
    double z1=0; 
    double z2=0;
    
    art::ServiceHandle<geo::Geometry> geo;
    
    geo->DetectorBigBox(&x1, &x2, &y1, &y2, &z1, &z2);
    
    if ( x1==0 && x2==0 && 
         y1==0 && y2==0 && 
         z1==0 && z2==0     ) {
      throw cet::exception("NoGeometryBoxCoords")
        << "No Geometry Box Coordinates Set\n"
        << __FILE__ << ":" << __LINE__ << "\n";
      return; 
    }
    
    // Loop over the particles stored by CRY
    int npart = truth.NParticles();
    for(int ipart = 0; ipart < npart; ++ipart){
      simb::MCParticle particle = truth.GetParticle(ipart);
      
      double px = particle.Px() / particle.P();
      double py = particle.Py() / particle.P();
      double pz = particle.Pz() / particle.P();

      double vx = particle.EndX();
      double vy = particle.EndY();
      double vz = particle.EndZ();

      // For the near detector, move all the cosmics from their position in
      // the window, straight down, to 1cm above the detector big box. This
      // ensures that a much larger fraction of them pass the IntersectTheBox
      // cut below. Subsequently the ProjectCosmicsToSurface call will
      // backtrack them up to the surface. The flux remains correct throughout
      // this scheme due to the homogeneity of the flux in space and time.
      // It does get the correlations wrong, but the chances of two correlated
      // events making it to the ND are remote anyway.
      if(geo->DetId() == novadaq::cnv::kNEARDET){
        vy = y2+1;
      }

      double xyz[3]  = { vx,  vy,  vz};
      double dxyz[3] = { px,  py,  pz};


      // If intersecting the DetectorBigBox, add particle
      //currently, using the methods in CosmicsGen code until the geometry methods can be validated
      if(this->isIntersectTheBox(xyz, dxyz, x1, x2, y1, y2, z1, z2)){
        simb::MCParticle part(particle.TrackId(),
                              particle.PdgCode(),
                              particle.Process(),
                              particle.Mother(),
                              particle.Mass(),
                              particle.StatusCode());
        part.AddTrajectoryPoint(TLorentzVector(vx, vy, vz, particle.T()),
                                particle.Momentum());

        mctruth_new.Add(part);
      }
    }// end of loop over particles

    // reassign the MCTruth
    truth = mctruth_new;
  }

  ////////////////////////////////////////////////////////////////////////
  // Method to project the cosmic particles back to the surface for the near detector.
  // CRYHelper is setup such that it projects CRY particles, which generated
  // at the surface up to the end of the World volume.
  // This is because ndos, for instance, is above the surface.
  // CRYHelper is standard (in the EventGeneratorBase), used by other experiments.
  // So, we don't want to modify it.
  // But we do project cosmics back to the surface for the near detector
  // (maybe should do the same for the far???)
  void CosmicsGen::ProjectCosmicsToSurface(simb::MCTruth& truth){

    art::ServiceHandle<geo::Geometry> geo;
    if(geo->DetId() != novadaq::cnv::kNEARDET) return;

    simb::MCTruth truth_after_surface_projection;

    double x1=0;
    double x2=0;
    double y1=0;
    double y2=0;
    double z1=0;
    double z2=0;

    // Get the World box
    geo->WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);

    if ( x1==0 && x2==0 && 
         y1==0 && y2==0 &&
         z1==0 && z2==0    ) {
      throw cet::exception("NoWorldBoxCoords")
        << "No World Box Coordinates Set\n"
        << __FILE__ << ":" << __LINE__ << "\n";
      return; 
    }

    // Loop over the particles
    int npart = truth.NParticles();
    for(int ipart = 0; ipart < npart; ++ipart){
      simb::MCParticle particle = truth.GetParticle(ipart);

      // make a new particle from this particle, where the first point 
      // is at the surface
      simb::MCParticle patsurface(particle.TrackId(),
                                  particle.PdgCode(),
                                  particle.Process(),
                                  particle.Mother(),
                                  particle.Mass(),
                                  particle.StatusCode());

      double px = particle.Px() / particle.P();
      double py = particle.Py() / particle.P();
      double pz = particle.Pz() / particle.P();

      double vx = particle.Vx();
      double vy = particle.Vy();
      double vz = particle.Vz();

      double xyz[3]  = { vx,  vy,  vz};
      double dxyz[3] = { px,  py,  pz};

      double xyzo[3];

      this->ProjectToSurface(xyz, dxyz, x1, x2, y1, y2, z1, z2, geo->SurfaceY(), xyzo);

      // reset vertex to be at the surface rather than the edge of the world
      TLorentzVector pos(xyzo[0], xyzo[1], xyzo[2], particle.T());
      TLorentzVector mom(particle.Px(), particle.Py(), particle.Pz(), particle.E());

      patsurface.AddTrajectoryPoint(pos, mom);

      truth_after_surface_projection.Add(patsurface);
    }// end of loop over particles

    // redefine truth
    truth = truth_after_surface_projection;
  }

  ////////////////////////////////////////////////////////////////////////
  // Method to project muons generated at the near detector to the DetectorBigBox.
  // This is to get rid of the muon propagation through the rock, since
  // the secondary particles would probably not reach the detector.
  // Also we account for the muon's energy loss.
  // After the muon reaches the DetectorBigBox, we start propagating it through Geant.
  void CosmicsGen::ProjectMuonsToDetectorBigBox(simb::MCTruth& truth){

    art::ServiceHandle<geo::Geometry> geo;

    // Projecting muons only if BigBox cut was used
    // Otherwise, the muon could miss the BigBox.
    // We just don't want to deal with these kind of events.
    if( !fIsBigBoxUsed || geo->DetId() != novadaq::cnv::kNEARDET) return;

    simb::MCTruth truth_after_bigbox_projection;

    double x1=0;
    double x2=0;
    double y1=0;
    double y2=0;
    double z1=0; 
    double z2=0;
    
    geo->DetectorBigBox(&x1, &x2, &y1, &y2, &z1, &z2);
    
    if ( x1==0 && x2==0 &&
         y1==0 && y2==0 &&
         z1==0 && z2==0    ) {
      throw cet::exception("NoGeometryBoxCoords")
        << "No Geometry Box Coordinates Set\n"
        << __FILE__ << ":" << __LINE__ << "\n";
      return; 
    }

    z2 += fProjectedZOffset;
     
    // Loop over the particles stored by CRY
    int npart = truth.NParticles();
    for(int ipart = 0; ipart < npart; ++ipart){

      simb::MCParticle particle = truth.GetParticle(ipart);

      //make a new particle from this particle, where the first point has the correct
      //momentum after the rock if it is a muon
      simb::MCParticle pafterrock(particle.TrackId(),
                                  particle.PdgCode(),
                                  particle.Process(),
                                  particle.Mother(),
                                  particle.Mass(),
                                  particle.StatusCode());

      // Project only mu+ and mu-
      if(abs(particle.PdgCode()) ==  13 ){

        double px = particle.Px() / particle.P();
        double py = particle.Py() / particle.P();
        double pz = particle.Pz() / particle.P();

        double vx = particle.Vx();
        double vy = particle.Vy();
        double vz = particle.Vz();

        double xyz[3]  = { vx,  vy,  vz};
        double dxyz[3] = { px,  py,  pz};

        double xyzo[3];
        double dlout=0;
        this->ProjectToBox(xyz, dxyz, x1, x2, y1, y2, z1, z2, xyzo, dlout);

        // In GeV
        double m    = particle.Mass();
        double etot = particle.E();
        double ke   = etot - m;

        double X = 2.5 * dlout; //g/sm^2 these magic numbers come from Gaisser, 
        double ksi = 250000.;   //g/cm^2 Cosmic Rays and Particle Physics equation (6.17)
        double epsilon = 500.0; //GeV

        double KE_after_rock = (ke + epsilon)*exp(-X/ksi) - epsilon;
        if(KE_after_rock < 0)continue;
        double E_after_rock = KE_after_rock + m;
                        
        double P_after_rock = sqrt(E_after_rock*E_after_rock - m*m);
        double Px = particle.Px()/particle.P() * P_after_rock;
        double Py = particle.Py()/particle.P() * P_after_rock;
        double Pz = particle.Pz()/particle.P() * P_after_rock;

        TLorentzVector pos(xyzo[0], xyzo[1], xyzo[2], particle.T());
        TLorentzVector mom(Px, Py, Pz, E_after_rock);

        pafterrock.AddTrajectoryPoint(pos, mom);


      }// end of if that particle is mu- or mu+
      else{
        // Don't waste time tracking non-muons through geant when they likely
        // won't reach the detector. As well as the waste of time, they'd break
        // the empty spill detection.
        continue;
      }

      // Add the particle into the new stack
      truth_after_bigbox_projection.Add(pafterrock);

    }// end of loop over particles
         
    // redefine the truth
    truth = truth_after_bigbox_projection;
  }

  /////////////////////////////////////////////////////
  /////////////////////////////////////////////////////
  //Method to separate MCTruth into separate objects for each particle
  void CosmicsGen::SeparateMCTruth(simb::MCTruth &truth,
                                   std::vector<simb::MCTruth>* truthcol){

    // Loop over the particles stored by CRY
    int npart = truth.NParticles();
    for(int ipart = 0; ipart < npart; ++ipart){
      simb::MCParticle particle = truth.GetParticle(ipart);

      simb::MCTruth ptruth;

      simb::MCParticle part(particle.TrackId(),
                            particle.PdgCode(),
                            particle.Process(),
                            particle.Mother(),
                            particle.Mass(),
                            particle.StatusCode());
      part.AddTrajectoryPoint(particle.Position(),
                              particle.Momentum());

      part.SetWeight(particle.Weight());

      ptruth.Add(part);

      truthcol->push_back(ptruth);
    }
   
  }

  //------------------------------------------------------------------------------
  // Suplementary methods
  ////////////////////////////////////////////////////////////////////////
  bool CosmicsGen::isIntersectTheBox(const double xyz[],
                                     const double dxyz[],
                                     double xlo, double xhi,
                                     double ylo, double yhi,
                                     double zlo, double zhi) {

    // If inside the box, obviously intesected
    if(xyz[0]>=xlo && xyz[0]<=xhi &&
       xyz[1]>=ylo && xyz[1]<=yhi &&
       xyz[2]>=zlo && xyz[2]<=zhi) return true;

    // So, now we know that the particle is outside the box

    double dx = 0., dy = 0., dz = 0.;

    // Checking intersection with 6 planes

    // 1. Check intersection with the upper plane
    dy = xyz[1] - yhi;
    // Check whether the track going from above and down or from below and up.
    // Otherwise not going to intersect this plane
    if( (dy > 0 && dxyz[1] < 0) ||
        (dy < 0 && dxyz[1] > 0) ){
      double dl = fabs( dy / dxyz[1] );
      double x = xyz[0] + dxyz[0] * dl;
      double z = xyz[2] + dxyz[2] * dl;

      // is it inside the square?
      if( x>=xlo && x<=xhi &&
          z>=zlo && z<=zhi) return true;
    }

    // 2. Check intersection with the lower plane
    dy = xyz[1] - ylo;
    // Check whether the track going from above and down or from below and up.
    // Otherwise not going to intersect this plane
    if( (dy > 0 && dxyz[1] < 0) ||
        (dy < 0 && dxyz[1] > 0) ){
      double dl = fabs( dy / dxyz[1] );
      double x = xyz[0] + dxyz[0] * dl;
      double z = xyz[2] + dxyz[2] * dl;

      // is it inside the square?
      if( x>=xlo && x<=xhi &&
          z>=zlo && z<=zhi) return true;
    }


    // 3. Check intersection with the right plane
    dz = xyz[2] - zhi;
    // Check whether the track going from right part to the left or from left part to right.
    // Otherwise not going to intersect this plane
    if( (dz > 0 && dxyz[2] < 0) ||
        (dz < 0 && dxyz[2] > 0) ){
      double dl = fabs( dz / dxyz[2] );
      double x = xyz[0] + dxyz[0] * dl;
      double y = xyz[1] + dxyz[1] * dl;

      // is it inside the square?
      if( x>=xlo && x<=xhi &&
          y>=ylo && y<=yhi) return true;
    }

    // 4. Check intersection with the left plane
    dz = xyz[2] - zlo;
    // Check whether the track going from right part to the left or from left part to right.
    // Otherwise not going to intersect this plane
    if( (dz > 0 && dxyz[2] < 0) ||
        (dz < 0 && dxyz[2] > 0) ){
      double dl = fabs( dz / dxyz[2] );
      double x = xyz[0] + dxyz[0] * dl;
      double y = xyz[1] + dxyz[1] * dl;

      // is it inside the square?
      if( x>=xlo && x<=xhi &&
          y>=ylo && y<=yhi) return true;
    }

    // 5. Check intersection with the far plane
    dx = xyz[0] - xhi;
    // Check whether the track going from far part toward us or from near part to right.
    // Otherwise not going to intersect this plane
    if( (dx > 0 && dxyz[0] < 0) ||
        (dx < 0 && dxyz[0] > 0) ){
      double dl = fabs( dx / dxyz[0] );
      double y = xyz[1] + dxyz[1] * dl;
      double z = xyz[2] + dxyz[2] * dl;

      // is it inside the square?
      if( z>=zlo && z<=zhi &&
          y>=ylo && y<=yhi) return true;
    }


    // 6. Check intersection with the near plane
    dx = xyz[0] - xlo;
    // Check whether the track going from far part toward us or from near part to right.
    // Otherwise not going to intersect this plane
    if( (dx > 0 && dxyz[0] < 0) ||
        (dx < 0 && dxyz[0] > 0) ){
      double dl = fabs( dx / dxyz[0] );
      double y = xyz[1] + dxyz[1] * dl;
      double z = xyz[2] + dxyz[2] * dl;

      // is it inside the square?
      if( z>=zlo && z<=zhi &&
          y>=ylo && y<=yhi) return true;
    }

    return false;
  }
  //
  ///----------------------------------------------------------------
  /// Project Back to the surface
  ///
  /// \param xyz        - The starting x,y,z location. Must be inside box.
  /// \param dxyz       - Direction vector
  /// \param xlo        - Low edge of box in x
  /// \param xhi        - High edge of box in x
  /// \param ylo        - Low edge of box in y
  /// \param yhi        - High edge of box in y
  /// \param zlo        - Low edge of box in z
  /// \param zhi        - High edge of box in z
  /// \param surface_y  - High edge of box in z
  /// \param xyzout     - On output, the position on the surface
  ///
  /// Note: It should be safe to use the same array for input and
  /// output.
  ///
  void CosmicsGen::ProjectToSurface(const double xyz[],
                                    const double dxyz[],
                                    double xlo, double xhi,
                                    double ylo, double yhi,
                                    double zlo, double zhi,
                                    const double surface_y,
                                    double xyzout[]) {
    art::ServiceHandle<geo::Geometry> geom;

    // Make sure that the cosmic muon is going down from above
    // (except in ND where we're actually rewinding them to the surface).
    assert(geom->DetId() == novadaq::cnv::kNEARDET ||
           xyz[1]  >= surface_y);
    assert(dxyz[1] <  0);
    // Also that the surface that we are projecting to is between ylo and yhi
    assert(surface_y >= ylo && surface_y <= yhi);

    // Step length
    double dl = - (xyz[1] - surface_y) / dxyz[1];

    // Make the step
    for (int i=0; i<3; ++i)
      xyzout[i] = xyz[i] + dxyz[i]*dl;


    if(geom->DetId() != novadaq::cnv::kNEARDET){
      // Make sure we're inside the box!
      assert(xyzout[0]>=xlo && xyzout[0]<=xhi);
      assert(xyzout[1]>=ylo && xyzout[1]<=yhi);
      assert(xyzout[2]>=zlo && xyzout[2]<=zhi);
    }
  }

  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  void CosmicsGen::ProjectToBox  (const double xyz[],
                                  const double dxyz[],
                                  double xlo, double xhi,
                                  double ylo, double yhi,
                                  double zlo, double zhi,
                                  double xyzout[], double& dlout) {


    // So, now we know that the particle is outside the box

    double dx, dy, dz;

    // Initialize to some big number
    dlout = 999999999.;

    // Checking intersection with 6 planes
    // 1. Check intersection with the upper plane
    dy = xyz[1] - yhi;
    // Check whether the track going from above and down or from below and up.
    // Otherwise not going to intersect this plane
    if( (dy > 0 && dxyz[1] < 0) ||
        (dy < 0 && dxyz[1] > 0) ){
      double dl = fabs( dy / dxyz[1] );
      double x = xyz[0] + dxyz[0] * dl;
      double z = xyz[2] + dxyz[2] * dl;

      // is it inside the square?
      if( x>=xlo && x<=xhi &&
          z>=zlo && z<=zhi &&
          dl < dlout ) {
        xyzout[0] = x;
        xyzout[1] = yhi;
        xyzout[2] = z;
        dlout = dl;
      }
    }

    // 2. Check intersection with the lower plane
    // Actually this cannot be the plane for the cosmics, because it would have to cross
    // some other plane first. But what the heck, let's check
    dy = xyz[1] - ylo;
    // Check whether the track going from above and down or from below and up.
    // Otherwise not going to intersect this plane
    if( (dy > 0 && dxyz[1] < 0) ||
        (dy < 0 && dxyz[1] > 0) ){
      double dl = fabs( dy / dxyz[1] );
      double x = xyz[0] + dxyz[0] * dl;
      double z = xyz[2] + dxyz[2] * dl;

      // is it inside the square?
      if( x>=xlo && x<=xhi &&
          z>=zlo && z<=zhi &&
          dl < dlout ) {
        xyzout[0] = x;
        xyzout[1] = ylo;
        xyzout[2] = z;
        dlout = dl;
      }
    }


    // 3. Check intersection with the right plane
    dz = xyz[2] - zhi;
    // Check whether the track going from right part to the left or from left part to right.
    // Otherwise not going to intersect this plane
    if( (dz > 0 && dxyz[2] < 0) ||
        (dz < 0 && dxyz[2] > 0) ){
      double dl = fabs( dz / dxyz[2] );
      double x = xyz[0] + dxyz[0] * dl;
      double y = xyz[1] + dxyz[1] * dl;

      // is it inside the square?
      if( x>=xlo && x<=xhi &&
          y>=ylo && y<=yhi &&
          dl < dlout ) {
        xyzout[0] = x;
        xyzout[1] = y;
        xyzout[2] = zhi;
        dlout = dl;
      }
    }

    // 4. Check intersection with the left plane
    dz = xyz[2] - zlo;
    // Check whether the track going from right part to the left or from left part to right.
    // Otherwise not going to intersect this plane
    if( (dz > 0 && dxyz[2] < 0) ||
        (dz < 0 && dxyz[2] > 0) ){
      double dl = fabs( dz / dxyz[2] );
      double x = xyz[0] + dxyz[0] * dl;
      double y = xyz[1] + dxyz[1] * dl;

      // is it inside the square?
      if( x>=xlo && x<=xhi &&
          y>=ylo && y<=yhi &&
          dl < dlout ) {
        xyzout[0] = x;
        xyzout[1] = y;
        xyzout[2] = zlo;
        dlout = dl;
      }
    }

    // 5. Check intersection with the far plane
    dx = xyz[0] - xhi;
    // Check whether the track going from far part toward us or from near part to right.
    // Otherwise not going to intersect this plane
    if( (dx > 0 && dxyz[0] < 0) ||
        (dx < 0 && dxyz[0] > 0) ){
      double dl = fabs( dx / dxyz[0] );
      double y = xyz[1] + dxyz[1] * dl;
      double z = xyz[2] + dxyz[2] * dl;

      // is it inside the square?
      if( z>=zlo && z<=zhi &&
          y>=ylo && y<=yhi &&
          dl < dlout ) {
        xyzout[0] = xhi;
        xyzout[1] = y;
        xyzout[2] = z;
        dlout = dl;
      }
    }


    // 6. Check intersection with the near plane
    dx = xyz[0] - xlo;
    // Check whether the track going from far part toward us or from near part to right.
    // Otherwise not going to intersect this plane
    if( (dx > 0 && dxyz[0] < 0) ||
        (dx < 0 && dxyz[0] > 0) ){
      double dl = fabs( dx / dxyz[0] );
      double y = xyz[1] + dxyz[1] * dl;
      double z = xyz[2] + dxyz[2] * dl;

      // is it inside the square?
      if( z>=zlo && z<=zhi &&
          y>=ylo && y<=yhi &&
          dl < dlout ) {
        xyzout[0] = xlo;
        xyzout[1] = y;
        xyzout[2] = z;
        dlout = dl;
      }
    }
  }

}// namespace

/////////////////////////////////////////////////////
namespace evgen
{
  DEFINE_ART_MODULE(CosmicsGen)
}