diff --git a/Analyses/src/CORSIKAGenPlots_module.cc b/Analyses/src/CORSIKAGenPlots_module.cc
new file mode 100644
index 0000000000..aa7fe8cd06
--- /dev/null
+++ b/Analyses/src/CORSIKAGenPlots_module.cc
@@ -0,0 +1,264 @@
+////////////////////////////////////////////////////////////////////////
+// Class:       CORSIKAGenPlots
+// Plugin Type: analyzer (art v2_10_04)
+// File:        CORSIKAGenPlots_module.cc
+//
+// Stefano Roberto Soleti (roberto@lbl.gov)
+//
+////////////////////////////////////////////////////////////////////////
+
+#include "art/Framework/Core/EDAnalyzer.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 "art_root_io/TFileService.h"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+
+#include "canvas/Utilities/InputTag.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+
+#include "GlobalConstantsService/inc/GlobalConstantsHandle.hh"
+#include "GlobalConstantsService/inc/PhysicsParams.hh"
+#include "GlobalConstantsService/inc/ParticleDataTable.hh"
+#include "MCDataProducts/inc/GenParticleCollection.hh"
+
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TTree.h"
+namespace mu2e {
+  class CORSIKAGenPlots;
+}
+
+using CLHEP::Hep3Vector;
+using CLHEP::HepLorentzVector;
+
+
+class mu2e::CORSIKAGenPlots : public art::EDAnalyzer {
+  public:
+    explicit CORSIKAGenPlots(fhicl::ParameterSet const & p);
+    // The compiler-generated destructor is fine for non-base
+    // classes without bare pointers or other resource use.
+
+    // Plugins should not be copied or assigned.
+    CORSIKAGenPlots(CORSIKAGenPlots const &) = delete;
+    CORSIKAGenPlots(CORSIKAGenPlots &&) = delete;
+    CORSIKAGenPlots & operator = (CORSIKAGenPlots const &) = delete;
+    CORSIKAGenPlots & operator = (CORSIKAGenPlots &&) = delete;
+
+    // Required functions.
+    void analyze(art::Event const & e) override;
+
+    // Selected optional functions.
+    void beginJob() override;
+    void endJob() override;
+
+  private:
+    std::string processName_;
+    std::string CORSIKAModuleLabel_;
+    std::string CORSIKAInstanceName_;
+    float _keMax = std::numeric_limits<float>::max();
+
+    // histograms
+    TH2F *_hXZ = nullptr;
+    TH1F *_hY = nullptr;
+    TH1F *_hKE = nullptr;
+    TH1F *_hTheta = nullptr;
+    TH1F *_hPhi = nullptr;
+    TH1F *_hPmag = nullptr;
+    TH1F *_hPyOverPmag = nullptr;
+    TH1F *_hTime = nullptr;
+    TH2F *_hPtypeKE = nullptr;
+    TH1F *_hNSecondaries = nullptr;
+
+    TTree *_cosmicTree = nullptr;
+
+    float _x = std::numeric_limits<float>::lowest();
+    float _y = std::numeric_limits<float>::lowest();
+    float _z = std::numeric_limits<float>::lowest();
+    float _px = std::numeric_limits<float>::lowest();
+    float _py = std::numeric_limits<float>::lowest();
+    float _pz = std::numeric_limits<float>::lowest();
+    float _theta = std::numeric_limits<float>::lowest();
+    float _phi = std::numeric_limits<float>::lowest();
+    float _KE = std::numeric_limits<float>::lowest();
+    float _p = std::numeric_limits<float>::lowest();
+    float _t = std::numeric_limits<float>::lowest();
+    PDGCode::type _pdgId;
+
+    void bookHists(art::ServiceHandle<art::TFileService> &);
+    GlobalConstantsHandle<ParticleDataTable> pdt;
+};
+
+mu2e::CORSIKAGenPlots::CORSIKAGenPlots(fhicl::ParameterSet const &p)
+    : EDAnalyzer(p),
+    processName_(p.get<std::string>("processName", "")),
+    CORSIKAModuleLabel_(p.get<std::string>("CORSIKAModuleLabel", "FromCorsikaBinary")),
+    CORSIKAInstanceName_(p.get<std::string>("CORSIKAInstanceName", "")),
+    _keMax(p.get<double>("keMax", 10E3))
+{
+  art::ServiceHandle<art::TFileService> tfs;
+
+  bookHists(tfs);
+
+  _cosmicTree = tfs->make<TTree>("cosmicTree", "TTree with cosmic ray info");
+
+  _cosmicTree->Branch("x", &_x, "x/F");
+  _cosmicTree->Branch("y", &_y, "y/F");
+  _cosmicTree->Branch("z", &_z, "z/F");
+  _cosmicTree->Branch("px", &_px, "px/F");
+  _cosmicTree->Branch("py", &_py, "py/F");
+  _cosmicTree->Branch("pz", &_pz, "pz/F");
+  _cosmicTree->Branch("theta", &_theta, "theta/F");
+  _cosmicTree->Branch("phi", &_phi, "phi/F");
+  _cosmicTree->Branch("KE", &_KE, "KE/F");
+  _cosmicTree->Branch("p", &_p, "p/F");
+  _cosmicTree->Branch("t", &_t, "t/F");
+  _cosmicTree->Branch("pdgId", &_pdgId, "pdgId/I");
+
+}
+
+
+void mu2e::CORSIKAGenPlots::analyze(art::Event const & e)
+{
+  art::Handle<GenParticleCollection> gpHandle;
+  bool success;
+
+  if (processName_.length() > 0)
+     success = e.getByLabel(CORSIKAModuleLabel_, CORSIKAInstanceName_, processName_,
+         gpHandle);
+  else
+    success = e.getByLabel(CORSIKAModuleLabel_, CORSIKAInstanceName_, gpHandle);
+
+  if (!success)
+    return;
+
+  const auto & particles = *gpHandle;
+
+
+  _hNSecondaries->Fill(particles.size());
+  for(const auto & p : particles)
+  {
+    _hXZ->Fill(p.position().x(), p.position().z());
+    _hY->Fill(p.position().y());
+
+    const HepPDT::ParticleData& p_data = pdt->particle(p.pdgId()).ref();
+    double mass = p_data.mass().value(); // in MeV
+
+    HepLorentzVector mom4 = p.momentum();
+    Hep3Vector mom3(mom4.px(), mom4.py(), mom4.pz());
+    Hep3Vector yMom3(mom4.pz(), mom4.px(), -mom4.py()); // for theta, phi
+
+    _hKE->Fill(mom4.e() - mass);
+    _hTheta->Fill(yMom3.theta());
+    _hPhi->Fill(yMom3.phi());
+
+    _hPmag->Fill(mom3.mag());
+    _hTime->Fill(p.time());
+    _hPyOverPmag->Fill(mom4.py() / mom3.mag());
+
+    _x = p.position().x();
+    _y = p.position().y();
+    _z = p.position().z();
+
+    _px = mom4.px();
+    _py = mom4.py();
+    _pz = mom4.pz();
+
+    _theta = yMom3.theta();
+    _phi = yMom3.phi();
+    _KE = mom4.e() - mass;
+    _p = mom3.mag();
+
+    _t = p.time();
+
+    _pdgId = p.pdgId();
+    _cosmicTree->Fill();
+
+    switch (p.pdgId()) {
+      case 13: // mu-
+        _hPtypeKE->Fill(mom4.e(), 0); break;
+      case -13: // mu+
+        _hPtypeKE->Fill(mom4.e(), 0); break;
+      case 22: // photon
+        _hPtypeKE->Fill(mom4.e(), 1); break;
+      case -11: // e+
+        _hPtypeKE->Fill(mom4.e(), 2); break;
+      case 11: // e-
+        _hPtypeKE->Fill(mom4.e(), 2); break;
+      case 2112: // neutron
+        _hPtypeKE->Fill(mom4.e(), 3); break;
+      case -2112: // neutron
+        _hPtypeKE->Fill(mom4.e(), 3); break;
+      case 2212: // proton
+        _hPtypeKE->Fill(mom4.e(), 4); break;
+      case -2212: // proton
+        _hPtypeKE->Fill(mom4.e(), 4); break;
+      case 111: // pi0
+        _hPtypeKE->Fill(mom4.e(), 5); break;
+      case 211: // pi+
+        _hPtypeKE->Fill(mom4.e(), 5); break;
+      case -211: // pi-
+        _hPtypeKE->Fill(mom4.e(), 5); break;
+      case 130: // k0 L
+        _hPtypeKE->Fill(mom4.e(), 6); break;
+      case 310: // k0 S
+        _hPtypeKE->Fill(mom4.e(), 6); break;
+      case 311: // k0
+        _hPtypeKE->Fill(mom4.e(), 6); break;
+      case 321: // k+
+        _hPtypeKE->Fill(mom4.e(), 6); break;
+      case -321: // k-
+        _hPtypeKE->Fill(mom4.e(), 6); break;
+      default: // others
+        _hPtypeKE->Fill(mom4.e(), 7); break;
+    }
+
+  }
+
+}
+
+void mu2e::CORSIKAGenPlots::beginJob()
+{
+  // Implementation of optional member function here.
+}
+
+void mu2e::CORSIKAGenPlots::endJob()
+{
+  // Implementation of optional member function here.
+}
+
+void mu2e::CORSIKAGenPlots::bookHists(art::ServiceHandle<art::TFileService> &tfs)
+{
+
+  _hXZ = tfs->make<TH2F>("XZ", "XZ", 500, -2.0e5,  2.0e5, 500, -2.0e5, 2.0e5 );
+  _hY = tfs->make<TH1F>("Y", "Y", 500, -15.0e3, 21.0e3 );
+
+  _hKE = tfs->make<TH1F>("E", "E", 2000, 0, _keMax);
+  _hTheta = tfs->make<TH1F>("Theta", "Theta",
+      200, -M_PI - 0.5, M_PI + 0.5);
+  _hPhi = tfs->make<TH1F>("Phi", "Phi",
+      200, -M_PI - 0.5, M_PI + 0.5);
+
+  _hPmag = tfs->make<TH1F>("Pmag", "Momentum modulus", 500, 0., 2E6);
+  _hPyOverPmag = tfs->make<TH1F>("PyOverPmag", "Py/Pmag", 200, -20., 20.);
+  _hTime = tfs->make<TH1F>("Time", "Timing of secondary particles",
+      1000, -1.E-3, 1E0);
+
+  _hNSecondaries = tfs->make<TH1F>("nSecondaries", "Number of secondaries",
+      1000, 0, 1000);
+
+  _hPtypeKE = tfs->make<TH2F>("PtypeKE", "Particle type vs energy",
+      2000, 0., 2E6, 8, 0, 8);
+  _hPtypeKE->GetYaxis()->SetBinLabel(1, "mu");
+  _hPtypeKE->GetYaxis()->SetBinLabel(2, "gamma");
+  _hPtypeKE->GetYaxis()->SetBinLabel(3, "electron");
+  _hPtypeKE->GetYaxis()->SetBinLabel(4, "neutron");
+  _hPtypeKE->GetYaxis()->SetBinLabel(5, "proton");
+  _hPtypeKE->GetYaxis()->SetBinLabel(6, "pion");
+  _hPtypeKE->GetYaxis()->SetBinLabel(7, "kaon");
+  _hPtypeKE->GetYaxis()->SetBinLabel(8, "others");
+}
+DEFINE_ART_MODULE(mu2e::CORSIKAGenPlots)
diff --git a/Analyses/src/CRYGenPlots_module.cc b/Analyses/src/CRYGenPlots_module.cc
index d2b65bdffd..d8ced9cd40 100644
--- a/Analyses/src/CRYGenPlots_module.cc
+++ b/Analyses/src/CRYGenPlots_module.cc
@@ -73,8 +73,23 @@ class mu2e::CRYGenPlots : public art::EDAnalyzer {
     TH2F *_hPtypeKE;
     TH1F *_hNSecondaries;
 
-    void bookHists();
-  GlobalConstantsHandle<ParticleDataTable> pdt;
+    TTree *_cosmicTree;
+
+    float _x;
+    float _y;
+    float _z;
+    float _px;
+    float _py;
+    float _pz;
+    float _theta;
+    float _phi;
+    float _KE;
+    float _p;
+    float _t;
+    int _pdgId;
+
+    void bookHists(art::ServiceHandle<art::TFileService> &);
+    GlobalConstantsHandle<ParticleDataTable> pdt;
 };
 
 
@@ -85,7 +100,24 @@ class mu2e::CRYGenPlots : public art::EDAnalyzer {
   , CRYInstanceName_(p.get<std::string>("CRYInstanceName", ""))
   , _keMax(p.get<double>("keMax", 10E3))
 {
-  bookHists();
+  art::ServiceHandle<art::TFileService> tfs;
+
+  bookHists(tfs);
+
+  _cosmicTree = tfs->make<TTree>("cosmicTree", "TTree with cosmic ray info");
+
+  _cosmicTree->Branch("x", &_x, "x/F");
+  _cosmicTree->Branch("y", &_y, "y/F");
+  _cosmicTree->Branch("z", &_z, "z/F");
+  _cosmicTree->Branch("px", &_px, "px/F");
+  _cosmicTree->Branch("py", &_py, "py/F");
+  _cosmicTree->Branch("pz", &_pz, "pz/F");
+  _cosmicTree->Branch("theta", &_theta, "theta/F");
+  _cosmicTree->Branch("phi", &_phi, "phi/F");
+  _cosmicTree->Branch("KE", &_KE, "KE/F");
+  _cosmicTree->Branch("p", &_p, "p/F");
+  _cosmicTree->Branch("t", &_t, "t/F");
+  _cosmicTree->Branch("pdgId", &_pdgId, "pdgId/I");
 }
 
 
@@ -115,7 +147,34 @@ void mu2e::CRYGenPlots::analyze(art::Event const & e)
 
     HepLorentzVector mom4 = p.momentum();
     Hep3Vector mom3(mom4.px(), mom4.py(), mom4.pz());
-    Hep3Vector yMom3(mom4.pz(), mom4.px(), mom4.py()); // for theta, phi
+    Hep3Vector yMom3(mom4.pz(), mom4.px(), -mom4.py()); // for theta, phi
+
+    _hKE->Fill(mom4.e() - mass);
+    _hTheta->Fill(yMom3.theta());
+    _hPhi->Fill(yMom3.phi());
+
+    _hPmag->Fill(mom3.mag());
+    _hTime->Fill(p.time());
+    _hPyOverPmag->Fill(mom4.py() / mom3.mag());
+
+    _x = p.position().x();
+    _y = p.position().y();
+    _z = p.position().z();
+
+    _px = mom4.px();
+    _py = mom4.py();
+    _pz = mom4.pz();
+
+    _theta = yMom3.theta();
+    _phi = yMom3.phi();
+
+    _KE = mom4.e() - mass;
+    _p = mom3.mag();
+
+    _t = p.time();
+
+    _pdgId = p.pdgId();
+    _cosmicTree->Fill();
 
     _hKE->Fill(mom4.e() - mass);
     _hTheta->Fill(yMom3.theta());
@@ -178,9 +237,8 @@ void mu2e::CRYGenPlots::endJob()
   // Implementation of optional member function here.
 }
 
-void mu2e::CRYGenPlots::bookHists()
+void mu2e::CRYGenPlots::bookHists(art::ServiceHandle<art::TFileService> &tfs)
 {
-  art::ServiceHandle<art::TFileService> tfs;
 
   _hXZ = tfs->make<TH2F>("XZ", "XZ", 500, -2.0e5,  2.0e5, 500, -2.0e5, 2.0e5 );
   _hY = tfs->make<TH1F>("Y", "Y", 500, -15.0e3, 21.0e3 );
diff --git a/EventGenerator/defaultConfigs/defaultCRYconfig.txt b/EventGenerator/defaultConfigs/defaultCRYconfig.txt
index 2ee87d2fd3..e753b6da50 100644
--- a/EventGenerator/defaultConfigs/defaultCRYconfig.txt
+++ b/EventGenerator/defaultConfigs/defaultCRYconfig.txt
@@ -36,7 +36,8 @@ int cosmicCRY.year = 2021;
 double cosmicCRY.latitude = 41.8;
 int cosmicCRY.altitude = 0; // meter, accepts either of 3 values: 0, 2100 or 11300
 double cosmicCRY.subboxLength = 300.; // meter
-double cosmicCRY.maxShowerEn  = 1e6; // meter
+double cosmicCRY.maxShowerEn  = 1e8; // MeV
+double cosmicCRY.minShowerEn  = 50; // MeV
 
 
 // This tells emacs to view this file in c++ mode.
diff --git a/EventGenerator/fcl/cryTest.fcl b/EventGenerator/fcl/cryTest.fcl
index 04ebb6d4ab..7c28ce3f0c 100644
--- a/EventGenerator/fcl/cryTest.fcl
+++ b/EventGenerator/fcl/cryTest.fcl
@@ -36,7 +36,8 @@ source : {
 services : {
   message : @local::default_message
   TFileService : { fileName : "nts.cryHist.owner.version.sequencer.root" }
-  GeometryService : { inputFile : "JobConfig/cosmic/geom_cosmic.txt" }
+  GeometryService : @local::Services.Core.GeometryService
+
   ConditionsService : { conditionsfile : "Mu2eG4/test/conditions_01.txt" }
   GlobalConstantsService : { inputFile : "Mu2eG4/test/globalConstants_01.txt" }
   G4Helper : { }
@@ -47,6 +48,9 @@ services : {
   }
 }
 
+services.GeometryService.inputFile : "JobConfig/cosmic/geom_cosmic.txt"
+
+
 physics : {
   producers : {
     cryGen : {
@@ -99,8 +103,8 @@ physics : {
     genCountLogger: { 
          module_type: GenEventCountReader 
     }
-    cryGenPlots: { 
-         module_type: CRYGenPlots 
+    cryGenPlots: {
+         module_type: CRYGenPlots
     }
   }
 
diff --git a/EventGenerator/inc/CosmicCRY.hh b/EventGenerator/inc/CosmicCRY.hh
index 2903319412..70bbd867b7 100644
--- a/EventGenerator/inc/CosmicCRY.hh
+++ b/EventGenerator/inc/CosmicCRY.hh
@@ -51,6 +51,7 @@ namespace mu2e {
       int _altitude;
       double _subboxLength;
       double _maxShowerEn;
+      double _minShowerEn;
 
       std::string _setupString;
       std::string _cryDataPath;
diff --git a/EventGenerator/src/CORSIKAEventGenerator_module.cc b/EventGenerator/src/CORSIKAEventGenerator_module.cc
new file mode 100644
index 0000000000..a3abcb58c6
--- /dev/null
+++ b/EventGenerator/src/CORSIKAEventGenerator_module.cc
@@ -0,0 +1,207 @@
+// Mu2e includes.
+#include "ConfigTools/inc/SimpleConfig.hh"
+#include "ConfigTools/inc/requireUniqueKey.hh"
+#include "MCDataProducts/inc/GenId.hh"
+#include "MCDataProducts/inc/GenParticleCollection.hh"
+#include "MCDataProducts/inc/G4BeamlineInfoCollection.hh"
+#include "SeedService/inc/SeedService.hh"
+
+// Includes from art and its toolchain.
+#include "art/Framework/Core/EDProducer.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Core/ModuleMacros.h"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+#include "art/Framework/Principal/Handle.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+
+#include "GlobalConstantsService/inc/GlobalConstantsHandle.hh"
+#include "GlobalConstantsService/inc/ParticleDataTable.hh"
+
+// C++ includes.
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <memory>
+
+#include "GeometryService/inc/GeomHandle.hh"
+#include "GeometryService/inc/GeometryService.hh"
+#include "GeometryService/inc/WorldG4.hh"
+#include "GeometryService/inc/DetectorSystem.hh"
+#include "GeometryService/inc/Mu2eEnvelope.hh"
+#include "MCDataProducts/inc/GenParticle.hh"
+#include "MCDataProducts/inc/GenParticleCollection.hh"
+#include "CalorimeterGeom/inc/Calorimeter.hh"
+#include "ExtinctionMonitorFNAL/Geometry/inc/ExtMonFNAL.hh"
+#include "GeneralUtilities/inc/safeSqrt.hh"
+
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "CLHEP/Vector/ThreeVector.h"
+#include "CLHEP/Random/RandFlat.h"
+
+#include "Mu2eUtilities/inc/VectorVolume.hh"
+
+
+using CLHEP::Hep3Vector;
+using CLHEP::HepLorentzVector;
+
+
+namespace mu2e {
+
+  class CorsikaEventGenerator : public art::EDProducer {
+    public:
+
+      struct Config
+      {
+        using Name = fhicl::Name;
+        using Comment = fhicl::Comment;
+
+        fhicl::Atom<std::string> corsikaModuleLabel{Name("corsikaModuleLabel"), Comment("Reference point in the coordinate system"), "FromCorsikaBinary"};
+        fhicl::Atom<std::string> refPointChoice{Name("refPointChoice"), Comment("Reference point in the coordinate system"), "UNDEFINED"};
+        fhicl::Atom<bool> projectToTargetBox{Name("projectToTargetBox"), Comment("Store only events that cross the target box"), false};
+        fhicl::Atom<float> targetBoxXmin{Name("targetBoxXmin"), Comment("Extension of the generation plane"), -5000};
+        fhicl::Atom<float> targetBoxXmax{Name("targetBoxXmax"), Comment("Extension of the generation plane"), 5000};
+        fhicl::Atom<float> targetBoxYmin{Name("targetBoxYmin"), Comment("Extension of the generation plane"), -5000};
+        fhicl::Atom<float> targetBoxYmax{Name("targetBoxYmax"), Comment("Extension of the generation plane"), 5000};
+        fhicl::Atom<float> targetBoxZmin{Name("targetBoxZmin"), Comment("Extension of the generation plane"), -5000};
+        fhicl::Atom<float> targetBoxZmax{Name("targetBoxZmax"), Comment("Extension of the generation plane"), 5000};
+      };
+      typedef art::EDProducer::Table<Config> Parameters;
+
+      explicit CorsikaEventGenerator(const Parameters &conf);
+      // Accept compiler written d'tor.  Modules are never moved or copied.
+      virtual void produce (art::Event& e);
+      virtual void beginSubRun(art::SubRun &sr);
+
+    private:
+
+      std::vector<CLHEP::Hep3Vector> _targetBoxIntersections;
+      std::vector<CLHEP::Hep3Vector> _worldIntersections;
+
+      bool  _geomInfoObtained = false;
+      Config _conf;
+      std::string _corsikaModuleLabel =  "FromCorsikaBinary";
+
+      float _worldXmin = 0;
+      float _worldXmax = 0;
+      float _worldYmin = 0;
+      float _worldYmax = 0;
+      float _worldZmin = 0;
+      float _worldZmax = 0;
+      float _targetBoxXmin = 0;
+      float _targetBoxXmax = 0;
+      float _targetBoxYmin = 0;
+      float _targetBoxYmax = 0;
+      float _targetBoxZmin = 0;
+      float _targetBoxZmax = 0;
+      bool _projectToTargetBox = false;
+      Hep3Vector _cosmicReferencePointInMu2e;
+      std::string _refPointChoice;
+
+  };
+
+  CorsikaEventGenerator::CorsikaEventGenerator(const Parameters &conf) : EDProducer{conf},
+                                                                         _conf(conf()),
+                                                                         _corsikaModuleLabel(_conf.corsikaModuleLabel()),
+                                                                         _targetBoxXmin(_conf.targetBoxXmin()),
+                                                                         _targetBoxXmax(_conf.targetBoxXmax()),
+                                                                         _targetBoxYmin(_conf.targetBoxYmin()),
+                                                                         _targetBoxYmax(_conf.targetBoxYmax()),
+                                                                         _targetBoxZmin(_conf.targetBoxZmin()),
+                                                                         _targetBoxZmax(_conf.targetBoxZmax()),
+                                                                         _projectToTargetBox(_conf.projectToTargetBox()),
+                                                                         _refPointChoice(_conf.refPointChoice())
+  {
+    produces<GenParticleCollection>();
+  }
+
+
+  void CorsikaEventGenerator::beginSubRun( art::SubRun &subrun){
+
+  }
+
+  void CorsikaEventGenerator::produce(art::Event &evt)
+  {
+    if (!_geomInfoObtained) {
+      GeomHandle<Mu2eEnvelope> env;
+      GeomHandle<WorldG4> worldGeom;
+      GeomHandle<DetectorSystem> detsys;
+
+      const float deltaX = 1; // mm
+      _worldXmin = worldGeom->mu2eOriginInWorld().x() - worldGeom->halfLengths()[0] + deltaX;
+      _worldXmax = worldGeom->mu2eOriginInWorld().x() + worldGeom->halfLengths()[0] - deltaX;
+      _worldYmin = env->ymin() + deltaX;
+      _worldYmax = env->ymax() - deltaX;
+      _worldZmin = worldGeom->mu2eOriginInWorld().z() - worldGeom->halfLengths()[2] + deltaX;
+      _worldZmax = worldGeom->mu2eOriginInWorld().z() + worldGeom->halfLengths()[2] - deltaX;
+
+      if (_refPointChoice == "TRACKER")
+      {
+        _cosmicReferencePointInMu2e = Hep3Vector(detsys->getOrigin().x(),
+                                                 0, detsys->getOrigin().z());
+      }
+      else if (_refPointChoice == "EXTMONFNAL")
+      {
+        GeomHandle<ExtMonFNAL::ExtMon> extMonFNAL;
+        _cosmicReferencePointInMu2e =
+            Hep3Vector(extMonFNAL->detectorCenterInMu2e().x(),
+                       0, extMonFNAL->detectorCenterInMu2e().z());
+      }
+      else if (_refPointChoice == "CALO")
+      {
+        GeomHandle<Calorimeter> calorimeter;
+        _cosmicReferencePointInMu2e = Hep3Vector(detsys->getOrigin().x(),
+                                                 0, calorimeter->disk(0).geomInfo().origin().z());
+      }
+      else if (_refPointChoice == "UNDEFINED")
+        _cosmicReferencePointInMu2e = Hep3Vector(0., 0, 0.);
+
+      _geomInfoObtained = true;
+    }
+
+    art::Handle<mu2e::GenParticleCollection> corsikaParticles;
+    evt.getByLabel(_corsikaModuleLabel, corsikaParticles);
+    const GenParticleCollection &particles(*corsikaParticles);
+    std::unique_ptr<GenParticleCollection> genParticles(new GenParticleCollection);
+
+    for (GenParticleCollection::const_iterator i = particles.begin(); i != particles.end(); ++i)
+    {
+      GenParticle const &particle = *i;
+      const HepLorentzVector mom4 = particle.momentum();
+
+      const Hep3Vector position(
+          particle.position().x() + _cosmicReferencePointInMu2e.x(),
+          particle.position().y(),
+          particle.position().z() + _cosmicReferencePointInMu2e.z());
+
+      if (_projectToTargetBox)
+      {
+        _worldIntersections.clear();
+        VectorVolume particleWorld(particle.position(), particle.momentum().vect(),
+                                   _worldXmin, _worldXmax, _worldYmin, _worldYmax, _worldZmin, _worldZmax);
+        particleWorld.calIntersections(_worldIntersections);
+
+        if (_worldIntersections.size() > 0)
+        {
+          // Being inside the world volume, the intersection can be only one.
+          const Hep3Vector projectedPos = _worldIntersections.at(0);
+          genParticles->push_back(GenParticle(static_cast<PDGCode::type>(particle.pdgId()),
+                                              GenId::cosmicCORSIKA, projectedPos, mom4,
+                                              particle.time()));
+        }
+      } else {
+        genParticles->push_back(GenParticle(static_cast<PDGCode::type>(particle.pdgId()),
+                  GenId::cosmicCORSIKA, position, mom4,
+                  particle.time()));
+      }
+    }
+
+    evt.put(std::move(genParticles));
+  }
+
+}
+
+
+using mu2e::CorsikaEventGenerator;
+DEFINE_ART_MODULE(CorsikaEventGenerator);
diff --git a/EventGenerator/src/CosmicCRY.cc b/EventGenerator/src/CosmicCRY.cc
index 0f746cd963..274e975988 100644
--- a/EventGenerator/src/CosmicCRY.cc
+++ b/EventGenerator/src/CosmicCRY.cc
@@ -42,369 +42,370 @@
 using CLHEP::Hep3Vector;
 using CLHEP::HepLorentzVector;
 
-namespace mu2e 
+namespace mu2e
 {
-  // The following part is needed for the RNG
-  template<class T> class RNGWrapper {
-    public:
-      static void set(T* object, double (T::*func)(void));
-      static double rng(void);
-    private:
-      static T* m_obj;
-      static double (T::*m_func)(void);
-  };// end of RNGWrapper class
-
-  template<class T> T* RNGWrapper<T>::m_obj;
-  template<class T> double (T::*RNGWrapper<T>::m_func)(void);
-  template<class T> void RNGWrapper<T>::set(T* object, double (T::*func)(void)) {
-    m_obj = object; m_func = func;
-  }
-  template<class T> double RNGWrapper<T>::rng(void) { return (m_obj->*m_func)(); }
-
-
-  CosmicCRY::CosmicCRY( art::Run& run,
-      const SimpleConfig& config, CLHEP::HepRandomEngine & engine )
-    : _verbose(config.getInt("cosmicCRY.verbose", 0) )
-        , _returnMuons( config.getBool("cosmicCRY.returnMuons", true))
-        , _returnNeutrons( config.getBool("cosmicCRY.returnNeutrons", true))
-        , _returnProtons( config.getBool("cosmicCRY.returnProtons", true))
-        , _returnGammas( config.getBool("cosmicCRY.returnGammas", true))
-        , _returnElectrons( config.getBool("cosmicCRY.returnElectrons", true))
-        , _returnPions( config.getBool("cosmicCRY.returnPions", true))
-        , _returnKaons( config.getBool("cosmicCRY.returnKaons", true))
-        , _month(config.getInt("cosmicCRY.month", 6))
-        , _day(config.getInt("cosmicCRY.day", 21))
-        , _year(config.getInt("cosmicCRY.year", 2020))
-        , _latitude( config.getDouble("cosmicCRY.latitude", 41.8))
-        , _altitude( config.getInt("cosmicCRY.altitude", 0))
-        , _subboxLength( config.getDouble("cosmicCRY.subboxLength", 100.))
-        , _maxShowerEn( config.getDouble("cosmicCRY.maxShowerEn", 1E6))
-        , _setupString("")
-        , _refY0(config.getDouble("cosmicCRY.refY0", 20000.))
-        , _refPointChoice(config.getString("cosmicCRY.refPoint", "UNDEFINED"))
-        , _directionChoice(config.getString("cosmicCRY.directionChoice", "ALL"))
-        , _cosmicReferencePointInMu2e()
-        , _vertical(false)
-        , _projectToTargetBox(config.getBool("cosmicCRY.projectToTargetBox", false))
-        , _geomInfoObtained(false)
-        , _targetBoxXmin(config.getDouble("cosmicCRY.targetBoxXmin", -1000))
-        , _targetBoxXmax(config.getDouble("cosmicCRY.targetBoxXmax", 1000))
-        , _targetBoxYmin(config.getDouble("cosmicCRY.targetBoxYmin", -1000))
-        , _targetBoxYmax(config.getDouble("cosmicCRY.targetBoxYmax", 1000))
-        , _targetBoxZmin(config.getDouble("cosmicCRY.targetBoxZmin", -1000))
-        , _targetBoxZmax(config.getDouble("cosmicCRY.targetBoxZmax", 1000))
-        {
-          _cryDataPath = std::string(std::getenv("CRYDATAPATH"));
-          if( _cryDataPath.length() == 0) 
-          {
-            mf::LogError("CosmicCRY") << "no variable CRYDATAPATH set. Exit now.";
-            throw cet::exception("Rethrow") << "This job cannot continue without a valid CRYDATAPATH.";
-          }
+// The following part is needed for the RNG
+template <class T>
+class RNGWrapper
+{
+public:
+  static void set(T *object, double (T::*func)(void));
+  static double rng(void);
+
+private:
+  static T *m_obj;
+  static double (T::*m_func)(void);
+}; // end of RNGWrapper class
+
+template <class T>
+T *RNGWrapper<T>::m_obj;
+template <class T>
+double (T::*RNGWrapper<T>::m_func)(void);
+template <class T>
+void RNGWrapper<T>::set(T *object, double (T::*func)(void))
+{
+  m_obj = object;
+  m_func = func;
+}
+template <class T>
+double RNGWrapper<T>::rng(void) { return (m_obj->*m_func)(); }
 
-          createSetupString();
-          _crySetup = new CRYSetup(_setupString, _cryDataPath);
+CosmicCRY::CosmicCRY(art::Run &run,
+                     const SimpleConfig &config, CLHEP::HepRandomEngine &engine)
+    : _verbose(config.getInt("cosmicCRY.verbose", 0)), _returnMuons(config.getBool("cosmicCRY.returnMuons", true)), _returnNeutrons(config.getBool("cosmicCRY.returnNeutrons", true)), _returnProtons(config.getBool("cosmicCRY.returnProtons", true)), _returnGammas(config.getBool("cosmicCRY.returnGammas", true)), _returnElectrons(config.getBool("cosmicCRY.returnElectrons", true)), _returnPions(config.getBool("cosmicCRY.returnPions", true)), _returnKaons(config.getBool("cosmicCRY.returnKaons", true)), _month(config.getInt("cosmicCRY.month", 6)), _day(config.getInt("cosmicCRY.day", 21)), _year(config.getInt("cosmicCRY.year", 2020)), _latitude(config.getDouble("cosmicCRY.latitude", 41.8)), _altitude(config.getInt("cosmicCRY.altitude", 0)), _subboxLength(config.getDouble("cosmicCRY.subboxLength", 100.)), _maxShowerEn(config.getDouble("cosmicCRY.maxShowerEn", 1E6)), _minShowerEn(config.getDouble("cosmicCRY.minShowerEn", 50)), _setupString(""), _refY0(config.getDouble("cosmicCRY.refY0", 20000.)), _refPointChoice(config.getString("cosmicCRY.refPoint", "UNDEFINED")), _directionChoice(config.getString("cosmicCRY.directionChoice", "ALL")), _cosmicReferencePointInMu2e(), _vertical(false), _projectToTargetBox(config.getBool("cosmicCRY.projectToTargetBox", false)), _geomInfoObtained(false), _targetBoxXmin(config.getDouble("cosmicCRY.targetBoxXmin", -1000)), _targetBoxXmax(config.getDouble("cosmicCRY.targetBoxXmax", 1000)), _targetBoxYmin(config.getDouble("cosmicCRY.targetBoxYmin", -1000)), _targetBoxYmax(config.getDouble("cosmicCRY.targetBoxYmax", 1000)), _targetBoxZmin(config.getDouble("cosmicCRY.targetBoxZmin", -1000)), _targetBoxZmax(config.getDouble("cosmicCRY.targetBoxZmax", 1000))
+{
+  _cryDataPath = std::string(std::getenv("CRYDATAPATH"));
+  if (_cryDataPath.length() == 0)
+  {
+    mf::LogError("CosmicCRY") << "no variable CRYDATAPATH set. Exit now.";
+    throw cet::exception("Rethrow") << "This job cannot continue without a valid CRYDATAPATH.";
+  }
 
-          RNGWrapper<CLHEP::HepRandomEngine>::set(
-              &engine,
-              &CLHEP::HepRandomEngine::flat);
-          _crySetup->setRandomFunction(RNGWrapper<CLHEP::HepRandomEngine>::rng);
+  createSetupString();
+  _crySetup = new CRYSetup(_setupString, _cryDataPath);
 
-          if (_verbose > 1) 
-            engine.showStatus();
+  RNGWrapper<CLHEP::HepRandomEngine>::set(
+      &engine,
+      &CLHEP::HepRandomEngine::flat);
+  _crySetup->setRandomFunction(RNGWrapper<CLHEP::HepRandomEngine>::rng);
 
-          _cryGen = std::make_unique<CRYGenerator>(_crySetup);
+  if (_verbose > 1)
+    engine.showStatus();
 
-          _GeV2MeV = CLHEP::GeV / CLHEP::MeV;
-          _m2mm = CLHEP::m / CLHEP::mm;
-	  _numEvents=0;
-        }
+  _cryGen = std::make_unique<CRYGenerator>(_crySetup);
 
+  _GeV2MeV = CLHEP::GeV / CLHEP::MeV;
+  _m2mm = CLHEP::m / CLHEP::mm;
+  _numEvents = 0;
+}
 
-  const double CosmicCRY::getLiveTime() {
-    return _cryGen->timeSimulated();
-  }
+const double CosmicCRY::getLiveTime()
+{
+  return _cryGen->timeSimulated();
+}
 
-  const double CosmicCRY::getShowerSumEnergy() {
-    return _showerSumEnergy;
-  }
+const double CosmicCRY::getShowerSumEnergy()
+{
+  return _showerSumEnergy;
+}
 
-  const unsigned long long int CosmicCRY::getNumEvents() {
-    return _numEvents;
-  }
+const unsigned long long int CosmicCRY::getNumEvents()
+{
+  return _numEvents;
+}
 
-  void CosmicCRY::generate( GenParticleCollection& genParts )
+void CosmicCRY::generate(GenParticleCollection &genParts)
+{
+  // Ref point, need to be here so that geometry info can be obtained
+  // , but let's do this only once
+  if (!_geomInfoObtained)
   {
-    // Ref point, need to be here so that geometry info can be obtained
-    // , but let's do this only once
-    if (!_geomInfoObtained) {
-      GeomHandle<Mu2eEnvelope> env;
-      GeomHandle<WorldG4>  worldGeom;
-      GeomHandle<ExtMonFNAL::ExtMon> extMonFNAL;
-      GeomHandle<DetectorSystem> detsys;
-
-      // slightly smaller box to avoid rounding error problem if any
-      double deltaX = 1; // mm
-      _worldXmin = worldGeom->mu2eOriginInWorld().x() - worldGeom->halfLengths()[0] + deltaX;
-      _worldXmax = worldGeom->mu2eOriginInWorld().x() + worldGeom->halfLengths()[0] - deltaX;
-      _worldYmin = env->ymin() + deltaX;
-      _worldYmax = env->ymax() - deltaX;
-      _worldZmin = worldGeom->mu2eOriginInWorld().z() - worldGeom->halfLengths()[2] + deltaX;
-      _worldZmax = worldGeom->mu2eOriginInWorld().z() + worldGeom->halfLengths()[2] - deltaX;
-
-      mf::LogInfo("CosmicCRY") << "Ref. point: " <<_cosmicReferencePointInMu2e << std::endl;
-      mf::LogInfo("CosmicCRY") << "Target box: " << _targetBoxXmin << ", " << 
-        _targetBoxXmax << ", " << _targetBoxYmin << ", " << _targetBoxYmax << 
-        ", " << _targetBoxZmin << ", " << _targetBoxZmax;
-      mf::LogInfo("CosmicCRY") << "World: " << _worldXmin << ", " << 
-        _worldXmax << ", " << _worldYmin << ", " << _worldYmax << ", " 
-        << _worldZmin << ", " << _worldZmax;
-
-      if (_refPointChoice == "TRACKER") 
-        _cosmicReferencePointInMu2e = Hep3Vector(detsys->getOrigin().x(),
-            _refY0, detsys->getOrigin().z());
-      else if (_refPointChoice == "EXTMONFNAL") 
-        _cosmicReferencePointInMu2e = 
+    GeomHandle<Mu2eEnvelope> env;
+    GeomHandle<WorldG4> worldGeom;
+    GeomHandle<ExtMonFNAL::ExtMon> extMonFNAL;
+    GeomHandle<DetectorSystem> detsys;
+
+    // slightly smaller box to avoid rounding error problem if any
+    double deltaX = 1; // mm
+    _worldXmin = worldGeom->mu2eOriginInWorld().x() - worldGeom->halfLengths()[0] + deltaX;
+    _worldXmax = worldGeom->mu2eOriginInWorld().x() + worldGeom->halfLengths()[0] - deltaX;
+    _worldYmin = env->ymin() + deltaX;
+    _worldYmax = env->ymax() - deltaX;
+    _worldZmin = worldGeom->mu2eOriginInWorld().z() - worldGeom->halfLengths()[2] + deltaX;
+    _worldZmax = worldGeom->mu2eOriginInWorld().z() + worldGeom->halfLengths()[2] - deltaX;
+
+    mf::LogInfo("CosmicCRY") << "Ref. point: " << _cosmicReferencePointInMu2e << std::endl;
+    mf::LogInfo("CosmicCRY") << "Target box: " << _targetBoxXmin << ", " << _targetBoxXmax << ", " << _targetBoxYmin << ", " << _targetBoxYmax << ", " << _targetBoxZmin << ", " << _targetBoxZmax;
+    mf::LogInfo("CosmicCRY") << "World: " << _worldXmin << ", " << _worldXmax << ", " << _worldYmin << ", " << _worldYmax << ", "
+                             << _worldZmin << ", " << _worldZmax;
+
+    if (_refPointChoice == "TRACKER")
+      _cosmicReferencePointInMu2e = Hep3Vector(detsys->getOrigin().x(),
+                                               _refY0, detsys->getOrigin().z());
+    else if (_refPointChoice == "EXTMONFNAL")
+      _cosmicReferencePointInMu2e =
           Hep3Vector(extMonFNAL->detectorCenterInMu2e().x(),
-              _refY0, extMonFNAL->detectorCenterInMu2e().z());
-      else if (_refPointChoice == "CALO") 
-      {
-        GeomHandle<Calorimeter> calorimeter;
-        _cosmicReferencePointInMu2e = Hep3Vector(detsys->getOrigin().x(),
-            _refY0, calorimeter->disk(0).geomInfo().origin().z());
-      }
-      else if (_refPointChoice == "UNDEFINED") 
-        _cosmicReferencePointInMu2e = Hep3Vector(0., _refY0, 0.);
-
-      _geomInfoObtained = true;
+                     _refY0, extMonFNAL->detectorCenterInMu2e().z());
+    else if (_refPointChoice == "CALO")
+    {
+      GeomHandle<Calorimeter> calorimeter;
+      _cosmicReferencePointInMu2e = Hep3Vector(detsys->getOrigin().x(),
+                                               _refY0, calorimeter->disk(0).geomInfo().origin().z());
     }
+    else if (_refPointChoice == "UNDEFINED")
+      _cosmicReferencePointInMu2e = Hep3Vector(0., _refY0, 0.);
 
-    bool passed = false;
-    // Getting CRY particles. Generate secondaries until you find one that intersects with the projected box
-    while(!passed){
-      std::vector<CRYParticle*> *secondaries = new std::vector<CRYParticle*>;
-      _cryGen->genEvent(secondaries);
-      _numEvents++;
-
-      double secondPtot = 0.;
-      _showerSumEnergy = 0.;
-
-      std::ostringstream oss;
-      for (unsigned j=0; j<secondaries->size(); j++) {
-	CRYParticle* secondary = (*secondaries)[j];
-
-	GlobalConstantsHandle<ParticleDataTable> pdt;
-	const HepPDT::ParticleData& p_data = pdt->particle(secondary->PDGid()).ref();
-	double mass = p_data.mass().value(); // in MeV
-
-
-	double ke = secondary->ke(); // MeV by default in CRY
-	double totalE = ke + mass;
-	_showerSumEnergy += totalE;
-	double totalP = safeSqrt(totalE * totalE - mass * mass);
-	
-	secondPtot += totalP;
-	
-	// Change coordinate system since y points upward, z points along
-	// the beam line; which make cry(xyz) -> mu2e(zxy), uvw -> mu2e(zxy)
-	Hep3Vector position(
-			    secondary->y() * _m2mm + _cosmicReferencePointInMu2e.x(),
-			    secondary->z() * _m2mm + _cosmicReferencePointInMu2e.y(),
-			    secondary->x() * _m2mm + _cosmicReferencePointInMu2e.z()); // to mm
-	HepLorentzVector mom4(totalP*secondary->v(), totalP*secondary->w(),
-			      totalP*secondary->u(), totalE);
-
-	if (_projectToTargetBox) {
-	  // Moving the CRY particle around: first find all intersections with
-	  // the target box, if there is any then find closest intersection with
-	  // world box
-	  _targetBoxIntersections.clear();
-	  calIntersections(position, mom4.vect(),
-			   _targetBoxIntersections, _targetBoxXmin, _targetBoxXmax,
-			   _targetBoxYmin, _targetBoxYmax, _targetBoxZmin, _targetBoxZmax);
-
-	  if (_targetBoxIntersections.size() > 0) {
-	  // if (_targetBoxIntersections.size() == 0) {
-	    _worldIntersections.clear();
-	    calIntersections(position, mom4.vect(), _worldIntersections,
-			     _worldXmin, _worldXmax, _worldYmin, _worldYmax, _worldZmin, _worldZmax);
-
-	    if (_worldIntersections.size() > 0) {
-	      int idx = 0;
-	      double closestDistance = distance(_worldIntersections.at(0), position);
-	      for (unsigned i = 0; i < _worldIntersections.size(); ++i) {
-		if (distance(_worldIntersections.at(i), position) < closestDistance) {
-		  idx = i;
-		  closestDistance = _targetBoxIntersections.at(idx).y();
-		}
-	      }
-	      
-	      Hep3Vector projectedPos = _worldIntersections.at(idx);
-	      // genParts.push_back(GenParticle(static_cast<PDGCode::type>(secondary->PDGid()),
-	      // 				     GenId::cosmicCRY, position, mom4,
-	      // 				     secondary->t() - _cryGen->timeSimulated()));
-	      genParts.push_back(GenParticle(static_cast<PDGCode::type>(secondary->PDGid()),
-	      				     GenId::cosmicCRY, projectedPos, mom4,
-	      				     secondary->t() - _cryGen->timeSimulated()));
-	    }
-	  }
-	}
-	else
-	  genParts.push_back(GenParticle(static_cast<PDGCode::type>(secondary->PDGid()),
-					 GenId::cosmicCRY, position, mom4,
-					 secondary->t() - _cryGen->timeSimulated()));
-
-	if (_verbose > 1) {
-	  oss << "Secondary " << j 
-	      << ": " << CRYUtils::partName(secondary->id()) 
-	      << " (pdgId " << secondary->PDGid() << ")"
-	      << ", kinetic energy " << secondary->ke()  << " MeV"
-	      << ", position " 
-	      << "(" << secondary->x()
-	      << ", " << secondary->y()
-	      << ", " << secondary->z()
-	      << ") m"
-	      << ", position at world boundaries" 
-	      << "("  << position.x()
-	      << ", " << position.y()
-	      << ", " << position.z()
-	      << ") m"
-	      << ", time " << secondary->t() << " sec"
-	      << ", mass: " << mass
-	      << ", mom: " << totalP
-	      << ", mom dir.: " << secondary->u() <<", " << secondary->v()
-	      << ", " << secondary->w() << "\n";
-	  
-	}
-	delete secondary;
-      }
-  
-      if (_verbose > 1) {
-	oss << "Total E: " << _showerSumEnergy;
-	mf::LogInfo("CosmicCRY") << oss.str();
-      }    
-      delete secondaries;
-
-      if(_showerSumEnergy<_maxShowerEn && genParts.size()>0)
-	passed = true;
-
-    }
+    _geomInfoObtained = true;
   }
 
-  void CosmicCRY::createSetupString()
+  bool passed = false;
+  // Getting CRY particles. Generate secondaries until you find one that intersects with the projected box
+  while (!passed)
   {
-    if (_returnMuons) 
-      _setupString.append("returnMuons 1 "); // must have the trailing white space
-    else
-      _setupString.append("returnMuons 0 ");
-
-    if (_returnNeutrons) 
-      _setupString.append("returnNeutrons 1 ");
-    else
-      _setupString.append("returnNeutrons 0 ");
-
-    if (_returnProtons) 
-      _setupString.append("returnProtons 1 ");
-    else
-      _setupString.append("returnProtons 0 ");
-
-    if (_returnGammas) 
-      _setupString.append("returnGammas 1 ");
-    else
-      _setupString.append("returnGammas 0 ");
-
-    if (_returnElectrons) 
-      _setupString.append("returnElectrons 1 ");
-    else
-      _setupString.append("returnElectrons 0 ");
-
-    if (_returnPions) 
-      _setupString.append("returnPions 1 ");
-    else
-      _setupString.append("returnPions 0 ");
-
-    if (_returnKaons) 
-      _setupString.append("returnKaons 1 ");
-    else
-      _setupString.append("returnKaons 0 ");
+    std::vector<CRYParticle *> *secondaries = new std::vector<CRYParticle *>;
+    _cryGen->genEvent(secondaries);
+    _numEvents++;
+
+    double secondPtot = 0.;
+    _showerSumEnergy = 0.;
 
     std::ostringstream oss;
-    oss << "date " << _month << "-" << _day << "-" << _year << " ";
-    oss << "latitude " << _latitude << " ";
-    oss << "altitude " << _altitude << " ";
-    oss << "subboxLength " << _subboxLength << " ";
+    for (unsigned j = 0; j < secondaries->size(); j++)
+    {
+      CRYParticle *secondary = (*secondaries)[j];
+
+      GlobalConstantsHandle<ParticleDataTable> pdt;
+      const HepPDT::ParticleData &p_data = pdt->particle(secondary->PDGid()).ref();
+      double mass = p_data.mass().value(); // in MeV
+
+      double ke = secondary->ke(); // MeV by default in CRY
+      if (ke < _minShowerEn)
+        continue;
+
+      double totalE = ke + mass;
+      _showerSumEnergy += totalE;
+      double totalP = safeSqrt(totalE * totalE - mass * mass);
+
+      secondPtot += totalP;
+
+      // Change coordinate system since y points upward, z points along
+      // the beam line; which make cry(xyz) -> mu2e(zxy), uvw -> mu2e(zxy)
+      Hep3Vector position(
+          secondary->y() * _m2mm + _cosmicReferencePointInMu2e.x(),
+          secondary->z() * _m2mm + _cosmicReferencePointInMu2e.y(),
+          secondary->x() * _m2mm + _cosmicReferencePointInMu2e.z()); // to mm
+      HepLorentzVector mom4(totalP * secondary->v(), totalP * secondary->w(),
+                            totalP * secondary->u(), totalE);
+
+      if (_projectToTargetBox)
+      {
+        // Moving the CRY particle around: first find all intersections with
+        // the target box, if there is any then find closest intersection with
+        // world box
+        _targetBoxIntersections.clear();
+        calIntersections(position, mom4.vect(),
+                         _targetBoxIntersections, _targetBoxXmin, _targetBoxXmax,
+                         _targetBoxYmin, _targetBoxYmax, _targetBoxZmin, _targetBoxZmax);
+
+        if (_targetBoxIntersections.size() > 0)
+        {
+          // if (_targetBoxIntersections.size() == 0) {
+          _worldIntersections.clear();
+          calIntersections(position, mom4.vect(), _worldIntersections,
+                           _worldXmin, _worldXmax, _worldYmin, _worldYmax, _worldZmin, _worldZmax);
 
-    _setupString.append(oss.str());
-  }
+          if (_worldIntersections.size() > 0)
+          {
+            int idx = 0;
+            double closestDistance = distance(_worldIntersections.at(0), position);
+            for (unsigned i = 0; i < _worldIntersections.size(); ++i)
+            {
+              if (distance(_worldIntersections.at(i), position) < closestDistance)
+              {
+                idx = i;
+                closestDistance = _targetBoxIntersections.at(idx).y();
+              }
+            }
+
+            Hep3Vector projectedPos = _worldIntersections.at(idx);
+            // genParts.push_back(GenParticle(static_cast<PDGCode::type>(secondary->PDGid()),
+            // 				     GenId::cosmicCRY, position, mom4,
+            // 				     secondary->t() - _cryGen->timeSimulated()));
+            genParts.push_back(GenParticle(static_cast<PDGCode::type>(secondary->PDGid()),
+                                           GenId::cosmicCRY, projectedPos, mom4,
+                                           secondary->t() - _cryGen->timeSimulated()));
+          }
+        }
+      }
+      else
+        genParts.push_back(GenParticle(static_cast<PDGCode::type>(secondary->PDGid()),
+                                       GenId::cosmicCRY, position, mom4,
+                                       secondary->t() - _cryGen->timeSimulated()));
 
-  void CosmicCRY::calIntersections(Hep3Vector orig, Hep3Vector dir,
-      std::vector<CLHEP::Hep3Vector> &intersections, double xMin, double xMax,
-      double yMin, double yMax, double zMin, double zMax) {
-    // roof: _targetBoxYmax, _targetBoxXmin, _targetBoxXmax, _targetBoxZmin, _targetBoxZmax
-    // skip projection if the particle goes parallely to the plane
-    if (dir.y() != 0.) {
-      double t = (yMax - orig.y()) / dir.y();
-      double x1 = dir.x() * t + orig.x();
-      double z1 = dir.z() * t + orig.z();
-      if (pointInBox(x1, z1, xMin, zMin, xMax, zMax)) {
-        intersections.push_back(Hep3Vector(x1, yMax, z1));
+      if (_verbose > 1)
+      {
+        oss << "Secondary " << j
+            << ": " << CRYUtils::partName(secondary->id())
+            << " (pdgId " << secondary->PDGid() << ")"
+            << ", kinetic energy " << secondary->ke() << " MeV"
+            << ", position "
+            << "(" << secondary->x()
+            << ", " << secondary->y()
+            << ", " << secondary->z()
+            << ") m"
+            << ", position at world boundaries"
+            << "(" << position.x()
+            << ", " << position.y()
+            << ", " << position.z()
+            << ") m"
+            << ", time " << secondary->t() << " sec"
+            << ", mass: " << mass
+            << ", mom: " << totalP
+            << ", mom dir.: " << secondary->u() << ", " << secondary->v()
+            << ", " << secondary->w() << "\n";
       }
+      delete secondary;
     }
 
-    //east: zMin, xMin, xMax, yMin, yMax
-    if (dir.z() != 0.) {
-      double t = (zMin - orig.z()) / dir.z();
-      double x1 = dir.x() * t + orig.x();
-      double y1 = dir.y() * t + orig.y();
-      if (pointInBox(x1, y1, xMin, yMin, xMax, yMax)) {
-        intersections.push_back(Hep3Vector(x1, y1, zMin));
-      }
+    if (_verbose > 1)
+    {
+      oss << "Total E: " << _showerSumEnergy;
+      mf::LogInfo("CosmicCRY") << oss.str();
     }
+    delete secondaries;
 
-    //west: zMax, xMin, xMax, yMin, yMax
-    if (dir.z() != 0.) {
-      double t = (zMax - orig.z()) / dir.z();
-      double x1 = dir.x() * t + orig.x();
-      double y1 = dir.y() * t + orig.y();
-      if (pointInBox(x1, y1, xMin, yMin, xMax, yMax)) {
-        intersections.push_back(Hep3Vector(x1, y1, zMax));
-      }
+    if (_showerSumEnergy < _maxShowerEn && genParts.size() > 0)
+      passed = true;
+  }
+}
+
+void CosmicCRY::createSetupString()
+{
+  if (_returnMuons)
+    _setupString.append("returnMuons 1 "); // must have the trailing white space
+  else
+    _setupString.append("returnMuons 0 ");
+
+  if (_returnNeutrons)
+    _setupString.append("returnNeutrons 1 ");
+  else
+    _setupString.append("returnNeutrons 0 ");
+
+  if (_returnProtons)
+    _setupString.append("returnProtons 1 ");
+  else
+    _setupString.append("returnProtons 0 ");
+
+  if (_returnGammas)
+    _setupString.append("returnGammas 1 ");
+  else
+    _setupString.append("returnGammas 0 ");
+
+  if (_returnElectrons)
+    _setupString.append("returnElectrons 1 ");
+  else
+    _setupString.append("returnElectrons 0 ");
+
+  if (_returnPions)
+    _setupString.append("returnPions 1 ");
+  else
+    _setupString.append("returnPions 0 ");
+
+  if (_returnKaons)
+    _setupString.append("returnKaons 1 ");
+  else
+    _setupString.append("returnKaons 0 ");
+
+  std::ostringstream oss;
+  oss << "date " << _month << "-" << _day << "-" << _year << " ";
+  oss << "latitude " << _latitude << " ";
+  oss << "altitude " << _altitude << " ";
+  oss << "subboxLength " << _subboxLength << " ";
+
+  _setupString.append(oss.str());
+}
+
+void CosmicCRY::calIntersections(Hep3Vector orig, Hep3Vector dir,
+                                 std::vector<CLHEP::Hep3Vector> &intersections, double xMin, double xMax,
+                                 double yMin, double yMax, double zMin, double zMax)
+{
+  // roof: _targetBoxYmax, _targetBoxXmin, _targetBoxXmax, _targetBoxZmin, _targetBoxZmax
+  // skip projection if the particle goes parallely to the plane
+  if (dir.y() != 0.)
+  {
+    double t = (yMax - orig.y()) / dir.y();
+    double x1 = dir.x() * t + orig.x();
+    double z1 = dir.z() * t + orig.z();
+    if (pointInBox(x1, z1, xMin, zMin, xMax, zMax))
+    {
+      intersections.push_back(Hep3Vector(x1, yMax, z1));
     }
+  }
 
-    //south: xMin, yMin, yMax, zMin, zMax
-    if (dir.x() != 0.) {
-      double t = (xMin - orig.x()) / dir.x();
-      double z1 = dir.z() * t + orig.z();
-      double y1 = dir.y() * t + orig.y();
-      if (pointInBox(z1, y1, zMin, yMin, zMax, yMax)) {
-        intersections.push_back(Hep3Vector(xMin, y1, z1));
-      }
+  //east: zMin, xMin, xMax, yMin, yMax
+  if (dir.z() != 0.)
+  {
+    double t = (zMin - orig.z()) / dir.z();
+    double x1 = dir.x() * t + orig.x();
+    double y1 = dir.y() * t + orig.y();
+    if (pointInBox(x1, y1, xMin, yMin, xMax, yMax))
+    {
+      intersections.push_back(Hep3Vector(x1, y1, zMin));
     }
+  }
 
-    //north: xMax, yMin, yMax, zMin, zMax
-    if (dir.x() != 0.) {
-      double t = (xMax - orig.x()) / dir.x();
-      double z1 = dir.z() * t + orig.z();
-      double y1 = dir.y() * t + orig.y();
-      if (pointInBox(z1, y1, zMin, yMin, zMax, yMax)) {
-        intersections.push_back(Hep3Vector(xMax, y1, z1));
-      }
+  //west: zMax, xMin, xMax, yMin, yMax
+  if (dir.z() != 0.)
+  {
+    double t = (zMax - orig.z()) / dir.z();
+    double x1 = dir.x() * t + orig.x();
+    double y1 = dir.y() * t + orig.y();
+    if (pointInBox(x1, y1, xMin, yMin, xMax, yMax))
+    {
+      intersections.push_back(Hep3Vector(x1, y1, zMax));
     }
+  }
 
+  //south: xMin, yMin, yMax, zMin, zMax
+  if (dir.x() != 0.)
+  {
+    double t = (xMin - orig.x()) / dir.x();
+    double z1 = dir.z() * t + orig.z();
+    double y1 = dir.y() * t + orig.y();
+    if (pointInBox(z1, y1, zMin, yMin, zMax, yMax))
+    {
+      intersections.push_back(Hep3Vector(xMin, y1, z1));
+    }
   }
 
-  bool CosmicCRY::pointInBox(double x, double y, double x0, double y0,
-      double x1, double y1)
+  //north: xMax, yMin, yMax, zMin, zMax
+  if (dir.x() != 0.)
   {
-    bool ret = false;
-    if ((x >= x0) && (x <= x1) && (y >= y0) && (y <= y1)) {
-      ret = true;
+    double t = (xMax - orig.x()) / dir.x();
+    double z1 = dir.z() * t + orig.z();
+    double y1 = dir.y() * t + orig.y();
+    if (pointInBox(z1, y1, zMin, yMin, zMax, yMax))
+    {
+      intersections.push_back(Hep3Vector(xMax, y1, z1));
     }
-    return ret;
   }
+}
 
-  double CosmicCRY::distance(const CLHEP::Hep3Vector &u, const CLHEP::Hep3Vector &v){
-    return safeSqrt((u.x() - v.x()) * (u.x() - v.x()) +
-        (u.y() - v.y()) * (u.y() - v.y()) +
-        (u.z() - v.z()) * (u.z() - v.z()));
+bool CosmicCRY::pointInBox(double x, double y, double x0, double y0,
+                           double x1, double y1)
+{
+  bool ret = false;
+  if ((x >= x0) && (x <= x1) && (y >= y0) && (y <= y1))
+  {
+    ret = true;
   }
+  return ret;
+}
+
+double CosmicCRY::distance(const CLHEP::Hep3Vector &u, const CLHEP::Hep3Vector &v)
+{
+  return safeSqrt((u.x() - v.x()) * (u.x() - v.x()) +
+                  (u.y() - v.y()) * (u.y() - v.y()) +
+                  (u.z() - v.z()) * (u.z() - v.z()));
 }
+} // namespace mu2e
\ No newline at end of file
diff --git a/EventGenerator/src/SConscript b/EventGenerator/src/SConscript
index 51e4768cb5..84d0f70314 100644
--- a/EventGenerator/src/SConscript
+++ b/EventGenerator/src/SConscript
@@ -17,6 +17,7 @@ crylib = os.environ['CRY_LIB']
 rootlibs  = env['ROOTLIBS']
 
 mainlib = helper.make_mainlib ( [
+        'mu2e_GlobalConstantsService',
         'mu2e_Mu2eUtilities',
         'mu2e_ConditionsService',
         'mu2e_GeometryService',
diff --git a/EventGenerator/test/corsikaTest.fcl b/EventGenerator/test/corsikaTest.fcl
new file mode 100644
index 0000000000..657c0fc6ce
--- /dev/null
+++ b/EventGenerator/test/corsikaTest.fcl
@@ -0,0 +1,147 @@
+#include "fcl/minimalMessageService.fcl"
+#include "fcl/standardProducers.fcl"
+#include "fcl/standardServices.fcl"
+
+
+process_name : corsikaTest
+
+targetParams: {
+    projectToTargetBox : true
+    targetBoxXmin: -10000
+    targetBoxXmax: 3000
+    targetBoxYmin: -5000
+    targetBoxYmax: 5000
+    targetBoxZmin: -5000
+    targetBoxZmax: 21000
+}
+
+source: {
+    module_type: FromCorsikaBinary
+    fileNames: ["/mu2e/app/users/srsoleti/corsika-77100/run/DAT110001"]
+    runNumber          : 1
+    firstSubRunNumber  : 0
+    firstEventNumber   : 1
+    showerAreaExtension  : 10000
+    maxEvents: -1
+    @table::targetParams
+    fluxConstant: 1.8e4
+}
+
+customCutPhoton: {
+   type: intersection
+   pars: [ { type: kineticEnergy cut: 10. },
+         { type: pdgId pars: [ 22 ] }]
+}
+customCutElectron: {
+   type: intersection
+   pars: [{type: kineticEnergy cut: 10. },
+         { type: pdgId pars: [ 11, -11 ] }]
+}
+customCutNeutron: {
+   type: intersection
+   // pars: [{type: kineticEnergy cut: 0.001},
+   pars: [{type: kineticEnergy cut: 10.000},
+         { type: pdgId pars: [ 2112, -2112 ] }]
+}
+
+
+services : {
+  message : @local::default_message
+  TFileService : { fileName : "nts.corsikaHist.owner.version.sequencer.root" }
+  GeometryService : @local::Services.Core.GeometryService
+
+  ConditionsService : { conditionsfile : "Mu2eG4/test/conditions_01.txt" }
+  GlobalConstantsService : { inputFile : "Mu2eG4/test/globalConstants_01.txt" }
+  G4Helper : { }
+  SeedService : @local::automaticSeeds
+  RandomNumberGenerator : { }
+
+  TimeTracker: {
+    printSummary: true  # default
+  }
+}
+
+services.GeometryService.inputFile : "JobConfig/cosmic/geom_cosmic.txt"
+
+physics : {
+  producers : {
+
+    corsikaGen : {
+      module_type : CORSIKAEventGenerator
+      corsikaModuleLabel: "FromCorsikaBinary"
+      refPointChoice: "UNDEFINED"
+      @table::targetParams
+    }
+
+    g4run :
+    {
+      module_type: Mu2eG4
+      physics: @local::mu2eg4DefaultPhysics
+      ResourceLimits: @local::mu2eg4DefaultResourceLimits
+	    TrajectoryControl: @local::mu2eg4DefaultTrajectories
+      debug: @local::mu2eg4DefaultDebug
+      visualization: @local::mu2eg4NoVisualization
+
+      generatorModuleLabel: corsikaGen
+      MultiStageParameters: {} // this is the first stage
+
+      SDConfig:
+      {
+         enableAllSDs: false
+         enableSD: [ CRV, calorimeter, stoppingtarget, tracker ]
+         TimeVD: { times: [] }
+      }
+
+      Mu2eG4StackingOnlyCut: @local::mu2eg4CutNeutrinos
+      Mu2eG4SteppingOnlyCut: {}
+      Mu2eG4CommonCut: {
+        type: union
+        pars: [@local::customCutPhoton,
+          @local::customCutElectron,
+          @local::customCutNeutron
+        ]
+        // write: customCut
+      }
+    }
+
+
+    randomsaver : @local::randomsaver
+
+    genCounter: {
+      module_type: GenEventCounter
+    }
+
+
+  }
+
+  analyzers : {
+
+    genCountLogger: {
+         module_type: GenEventCountReader 
+    }
+
+    corsikaGenPlots: {
+         CORSIKAModuleLabel   : "corsikaGen"
+         CORSIKAInstanceName  : ""
+         module_type: CORSIKAGenPlots
+    }
+
+  }
+
+  trigFilter: [corsikaGen, g4run, g4run, randomsaver]
+  e1 : [outfile]
+  ana : [corsikaGenPlots]
+
+  trigger_paths  : [trigFilter]
+  end_paths : [e1, ana]
+}
+
+outputs : {
+  outfile : {
+    module_type : RootOutput
+    fileName : "sim.corsikaOut.owner.version.sequencer.art"
+  }
+}
+services.SeedService.baseSeed         :  3425
+services.SeedService.maxUniqueEngines :  20
+// vim: set ft=cpp:
diff --git a/MCDataProducts/inc/GenId.hh b/MCDataProducts/inc/GenId.hh
index 6ef30da93c..999cb4af10 100644
--- a/MCDataProducts/inc/GenId.hh
+++ b/MCDataProducts/inc/GenId.hh
@@ -44,8 +44,8 @@ namespace mu2e {
       MARS, StoppedParticleReactionGun, bremElectronGun, muonicXRayGun, //30
       fromSimParticleStartPoint, fromSimParticleCompact, StoppedParticleG4Gun, //33
       CaloCalib, InFlightParticleSampler, muplusDecayGun, StoppedMuonXRayGammaRayGun, //37
-      cosmicCRY, pbarFlat, fromAscii, ExternalRMC, InternalRMC, CeLeadingLog, //43
-      lastEnum //44
+      cosmicCRY,  pbarFlat, fromAscii, ExternalRMC, InternalRMC, CeLeadingLog, cosmicCORSIKA, //44
+      lastEnum //45
     };
 
     // Keep this in sync with the enum. Used in GenId.cc
@@ -60,7 +60,7 @@ namespace mu2e {
       "MARS", "StoppedParticleReactionGun","bremElectronGun", "muonicXRayGun", \
       "fromSimParticleStartPoint", "fromSimParticleCompact", "StoppedParticleG4Gun", \
       "CaloCalib", "InFlightParticleSampler","muplusDecayGun", "StoppedMuonXRayGammaRayGun", \
-      "CosmicCRY","pbarFlat","fromAscii","ExternalRMC","InternalRMC","CeLeadingLog"
+      "CosmicCRY", "pbarFlat","fromAscii","ExternalRMC","InternalRMC","CeLeadingLog", "CosmicCORSIKA",
 
   public:
 
diff --git a/Mu2eUtilities/inc/VectorVolume.hh b/Mu2eUtilities/inc/VectorVolume.hh
new file mode 100644
index 0000000000..3b25a7f634
--- /dev/null
+++ b/Mu2eUtilities/inc/VectorVolume.hh
@@ -0,0 +1,52 @@
+#ifndef Mu2eUtilities_VectorVolume_hh
+#define Mu2eUtilities_VectorVolume_hh
+//
+// Given a vector and a volume, gives the intersections between the two
+//
+// Original author Stefano Roberto Soleti
+//
+//
+// Constructor arguments:
+
+//
+//
+
+// CLHEP includes
+#include "CLHEP/Vector/ThreeVector.h"
+#include <vector>
+#include "GeneralUtilities/inc/safeSqrt.hh"
+using namespace std;
+
+namespace mu2e {
+
+  class VectorVolume{
+
+  public:
+    VectorVolume(CLHEP::Hep3Vector const& _orig,
+                 CLHEP::Hep3Vector const& _dir,
+                 float _xMin, float _xMax,
+                 float _yMin, float _yMax,
+                 float _zMin, float _zMax);
+    ~VectorVolume();
+
+    bool pointInBox(float x, float y,
+                    float x0, float y0,
+                    float x1, float y1);
+    float distance(const CLHEP::Hep3Vector &u, const CLHEP::Hep3Vector &v);
+    void calIntersections(std::vector<CLHEP::Hep3Vector> & intersections);
+
+  private:
+    CLHEP::Hep3Vector _orig;
+    CLHEP::Hep3Vector _dir;
+    float _xMin;
+    float _xMax;
+    float _yMin;
+    float _yMax;
+    float _zMin;
+    float _zMax;
+
+  };
+
+} // namespace mu2e
+
+#endif /* Mu2eUtilities_VectorVolume_hh */
diff --git a/Mu2eUtilities/src/VectorVolume.cc b/Mu2eUtilities/src/VectorVolume.cc
new file mode 100644
index 0000000000..fbcdcc6b97
--- /dev/null
+++ b/Mu2eUtilities/src/VectorVolume.cc
@@ -0,0 +1,119 @@
+//
+// Given a vector and a volume, gives the intersections between the two
+//
+//
+// _original author Stefano Roberto Soleti
+//
+
+#include <iostream>
+#include "Mu2eUtilities/inc/VectorVolume.hh"
+
+using namespace std;
+
+namespace mu2e {
+
+
+  VectorVolume::VectorVolume(CLHEP::Hep3Vector const& orig,
+                             CLHEP::Hep3Vector const& dir,
+                             float xMin, float xMax,
+                             float yMin, float yMax,
+                             float zMin, float zMax):
+                             _orig(orig),
+                             _dir(dir),
+                             _xMin(xMin),
+                             _xMax(xMax),
+                             _yMin(yMin),
+                             _yMax(yMax),
+                             _zMin(zMin),
+                             _zMax(zMax)
+  {
+  }
+
+  VectorVolume::~VectorVolume(){
+  }
+
+  bool VectorVolume::pointInBox(float x, float y, float x0, float y0,
+                                        float x1, float y1)
+  {
+    bool ret = false;
+    if ((x >= x0) && (x <= x1) && (y >= y0) && (y <= y1))
+    {
+      ret = true;
+    }
+    return ret;
+  }
+
+  float VectorVolume::distance(const CLHEP::Hep3Vector &u, const CLHEP::Hep3Vector &v)
+  {
+    return safeSqrt((u.x() - v.x()) * (u.x() - v.x()) +
+                    (u.y() - v.y()) * (u.y() - v.y()) +
+                    (u.z() - v.z()) * (u.z() - v.z()));
+  }
+
+  void VectorVolume::calIntersections(std::vector<CLHEP::Hep3Vector> & intersections)
+  {
+    // roof: _targetBox_yMax, _targetBox_xMin, _targetBox_xMax, _targetBox_zMin, _targetBox_zMax
+    // skip projection if the particle goes parallely to the plane
+    if (_dir.y() != 0.)
+    {
+      const float t = (_yMax - _orig.y()) / _dir.y();
+      const float x1 = _dir.x() * t + _orig.x();
+      const float z1 = _dir.z() * t + _orig.z();
+      // std::cout << "x1 " << x1 << " " << z1 << std::endl;
+
+      if (pointInBox(x1, z1, _xMin, _zMin, _xMax, _zMax))
+      {
+        intersections.push_back(CLHEP::Hep3Vector(x1, _yMax, z1));
+      }
+    }
+
+    //east: _zMin, _xMin, _xMax, _yMin, _yMax
+    if (_dir.z() != 0.)
+    {
+      const float t = (_zMin - _orig.z()) / _dir.z();
+      const float x1 = _dir.x() * t + _orig.x();
+      const float y1 = _dir.y() * t + _orig.y();
+      if (pointInBox(x1, y1, _xMin, _yMin, _xMax, _yMax))
+      {
+        intersections.push_back(CLHEP::Hep3Vector(x1, y1, _zMin));
+      }
+    }
+
+    //west: _zMax, _xMin, _xMax, _yMin, _yMax
+    if (_dir.z() != 0.)
+    {
+      const float t = (_zMax - _orig.z()) / _dir.z();
+      const float x1 = _dir.x() * t + _orig.x();
+      const float y1 = _dir.y() * t + _orig.y();
+      if (pointInBox(x1, y1, _xMin, _yMin, _xMax, _yMax))
+      {
+        intersections.push_back(CLHEP::Hep3Vector(x1, y1, _zMax));
+      }
+    }
+
+    //south: _xMin, _yMin, _yMax, _zMin, _zMax
+    if (_dir.x() != 0.)
+    {
+      const float t = (_xMin - _orig.x()) / _dir.x();
+      const float z1 = _dir.z() * t + _orig.z();
+      const float y1 = _dir.y() * t + _orig.y();
+      if (pointInBox(z1, y1, _zMin, _yMin, _zMax, _yMax))
+      {
+        intersections.push_back(CLHEP::Hep3Vector(_xMin, y1, z1));
+      }
+    }
+
+    //north: _xMax, _yMin, _yMax, _zMin, _zMax
+    if (_dir.x() != 0.)
+    {
+      const float t = (_xMax - _orig.x()) / _dir.x();
+      const float z1 = _dir.z() * t + _orig.z();
+      const float y1 = _dir.y() * t + _orig.y();
+      if (pointInBox(z1, y1, _zMin, _yMin, _zMax, _yMax))
+      {
+        intersections.push_back(CLHEP::Hep3Vector(_xMax, y1, z1));
+      }
+    }
+  }
+
+} // end namespace mu2e
diff --git a/Sources/inc/CosmicCORSIKA.hh b/Sources/inc/CosmicCORSIKA.hh
new file mode 100644
index 0000000000..da94f6028a
--- /dev/null
+++ b/Sources/inc/CosmicCORSIKA.hh
@@ -0,0 +1,246 @@
+// Cosmic rays generator using CORSIKA
+
+#ifndef Sources_inc_CosmicCORSIKA_hh
+#define Sources_inc_CosmicCORSIKA_hh
+
+#include <vector>
+
+
+#include "GlobalConstantsService/inc/GlobalConstantsHandle.hh"
+#include "GlobalConstantsService/inc/ParticleDataTable.hh"
+
+// Framework includes
+#include "art/Framework/Core/EDProducer.h"
+#include "art/Framework/Core/ModuleMacros.h"
+#include "art/Framework/Principal/Event.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+
+#include "GeometryService/inc/GeomHandle.hh"
+#include "GeometryService/inc/GeometryService.hh"
+#include "GeometryService/inc/WorldG4.hh"
+#include "GeometryService/inc/DetectorSystem.hh"
+#include "GeometryService/inc/Mu2eEnvelope.hh"
+#include "MCDataProducts/inc/GenParticle.hh"
+#include "MCDataProducts/inc/GenParticleCollection.hh"
+#include "CalorimeterGeom/inc/Calorimeter.hh"
+#include "ExtinctionMonitorFNAL/Geometry/inc/ExtMonFNAL.hh"
+#include "GeneralUtilities/inc/safeSqrt.hh"
+
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "CLHEP/Vector/ThreeVector.h"
+#include "CLHEP/Random/RandFlat.h"
+
+#include "fhiclcpp/types/ConfigurationTable.h"
+
+#include "Mu2eUtilities/inc/VectorVolume.hh"
+
+
+namespace art
+{
+  class Run;
+}
+
+struct Config
+{
+  using Name = fhicl::Name;
+  using Comment = fhicl::Comment;
+  fhicl::Sequence<std::string> showerInputFiles{Name("fileNames"),Comment("List of CORSIKA binary output paths")};
+  fhicl::Atom<int> firstEventNumber{Name("firstEventNumber"), Comment("Reference point on the Y axis"), 0};
+  fhicl::Atom<int> firstSubRunNumber{Name("firstSubRunNumber"), Comment("Reference point on the Y axis"), 0};
+  fhicl::Atom<int> maxEvents{Name("maxEvents"), Comment("Reference point on the Y axis"), 0};
+  fhicl::Atom<unsigned int> runNumber{Name("runNumber"), Comment("Reference point on the Y axis"), 0};
+  fhicl::Atom<std::string> module_label{Name("module_label"), Comment("Reference point on the Y axis"), ""};
+  fhicl::Atom<std::string> module_type{Name("module_type"), Comment("Reference point on the Y axis"), ""};
+  fhicl::Atom<bool> projectToTargetBox{Name("projectToTargetBox"), Comment("Store only events that cross the target box"), false};
+  fhicl::Atom<float> showerAreaExtension{Name("showerAreaExtension"), Comment("Reference point on the Y axis"), 10000};
+  fhicl::Atom<float> tOffset{Name("tOffset"), Comment("Time offset"), 0};
+  fhicl::Atom<float> fluxConstant{Name("fluxConstant"), Comment("Flux constant"), 1.8e4};
+  fhicl::Atom<float> targetBoxXmin{Name("targetBoxXmin"), Comment("Extension of the generation plane"), -5000};
+  fhicl::Atom<float> targetBoxXmax{Name("targetBoxXmax"), Comment("Extension of the generation plane"), 5000};
+  fhicl::Atom<float> targetBoxYmin{Name("targetBoxYmin"), Comment("Extension of the generation plane"), -5000};
+  fhicl::Atom<float> targetBoxYmax{Name("targetBoxYmax"), Comment("Extension of the generation plane"), 5000};
+  fhicl::Atom<float> targetBoxZmin{Name("targetBoxZmin"), Comment("Extension of the generation plane"), -5000};
+  fhicl::Atom<float> targetBoxZmax{Name("targetBoxZmax"), Comment("Extension of the generation plane"), 5000};
+
+};
+typedef fhicl::WrappedTable<Config> Parameters;
+
+class GenParticleCollection;
+
+namespace mu2e {
+
+  class CosmicCORSIKA{
+
+    public:
+      CosmicCORSIKA(const Config& conf);
+      // CosmicCORSIKA(art::Run &run, CLHEP::HepRandomEngine &engine);
+      ~CosmicCORSIKA();
+      const float getLiveTime();
+      const unsigned int getNumShowers();
+
+      const std::map<unsigned int, int> corsikaToPdgId = {
+        {1, 22}, // gamma
+        {2, -11}, // e+
+        {3, 11}, // e-
+        {5, -13}, // mu+
+        {6, 13}, // mu-
+        {7, 111}, // pi0
+        {8, 211}, // pi+
+        {9, -211}, // pi-
+        {10, 130}, // K0_L
+        {11, 321}, // K+
+        {12, -321}, // K-
+        {13, 2112}, // n
+        {14, 2212}, // p
+        {15, -2212}, // pbar
+        {16, 310}, // K0_S
+        {17, 221}, // eta
+        {18, 3122}, // Lambda
+        {19, 3222}, // Sigma+
+        {20, 3212}, // Sigma0
+        {21, 3112}, // Sigma-
+        {22, 3322}, // Cascade0
+        {23, 3312}, // Cascade-
+        {24, 3334}, // Omega-
+        {25, -2112}, // nbar
+        {26, -3122}, // Lambdabar
+        {27, -3112}, // Sigma-bar
+        {28, -3212}, // Sigma0bar
+        {29, -3222}, // Sigma+bar
+        {30, -3322}, // Cascade0bar
+        {31, -3312}, // Cascade+bar
+        {32, -3334}, // Omega+bar
+
+        {50, 223}, // omega
+        {51, 113}, // rho0
+        {52, 213}, // rho+
+        {53, -213}, // rho-
+        {54, 2224}, // Delta++
+        {55, 2214}, // Delta+
+        {56, 2114}, // Delta0
+        {57, 1114}, // Delta-
+        {58, -2224}, // Delta--bar
+        {59, -2214}, // Delta-bar
+        {60, -2114}, // Delta0bar
+        {61, -1114}, // Delta+bar
+        {62, 10311}, // K*0
+        {63, 10321}, // K*+
+        {64, -10321}, // K*-
+        {65, -10311}, // K*0bar
+        {66, 12}, // nu_e
+        {67, -12}, // nu_ebar
+        {68, 14}, // nu_mu
+        {69, -14}, // nu_mubar
+
+        {116, 421}, // D0
+        {117, 411}, // D+
+        {118, -411}, // D-bar
+        {119, -421}, // D0bar
+        {120, 431}, // D+_s
+        {121, -431}, // D-_sbar
+        {122, 441}, // eta_c
+        {123, 423}, // D*0
+        {124, 413}, // D*+
+        {125, -413}, // D*-bar
+        {126, -423}, // D*0bar
+        {127, 433}, // D*+_s
+        {128, -433}, // D*-_s
+
+        {130, 443}, // J/Psi
+        {131, -15}, // tau+
+        {132, 15}, // tau-
+        {133, 16}, // nu_tau
+        {134, -16}, // nu_taubar
+
+        {137, 4122}, // Lambda+_c
+        {138, 4232}, // Cascade+_c
+        {139, 4132}, // Cascade0_c
+        {140, 4222}, // Sigma++_c
+        {141, 4212}, // Sigma+_c
+        {142, 4112}, // Sigma0_c
+        {143, 4322}, // Cascade'+_c
+        {144, 4312}, // Cascade'0_c
+        {145, 4332}, // Omega0_c
+        {149, -4122}, // Lambda-_cbar
+        {150, -4232}, // Cascade-_cbar
+        {151, -4132}, // Cascade0_cbar
+        {152, -4222}, // Sigma--_cbar
+        {153, -4212}, // Sigma-_cbar
+        {154, -4112}, // Sigma0_cbar
+        {155, -4322}, // Cascade'-_cbar
+        {156, -4312}, // Cascade'0_cbar
+        {157, -4332}, // Omega0_cbar
+        {161, 4224}, // Sigma*++_c
+        {162, 1214}, // Sigma*+_c
+        {163, 4114}, // Sigma*0_c
+
+        {171, -4224}, // Sigma*--_cbar
+        {172, -1214}, // Sigma*-_cbar
+        {173, -4114}, // Sigma*0_cbar
+        {176, 511}, // B0
+        {177, 521}, // B+
+        {178, -521}, // B-bar
+        {179, -511}, // B0bar
+        {180, 531}, // B0_s
+        {181, -531}, // B0_sbar
+        {182, 541}, // B+_c
+        {183, -541}, // B-_cbar
+        {184, 5122}, // Lambda0_b
+        {185, 5112}, // Sigma-_b
+        {186, 5222}, // Sigma+_b
+        {187, 5232}, // Cascade0_b
+        {188, 5132}, // Cascade-_b
+        {189, 5332}, // Omega-_b
+        {190, -5112}, // Lambda0_bbar
+        {191, -5222}, // Sigma+_bbar
+        {192, -5112}, // Sigma-_bbar
+        {193, -5232}, // Cascade0_bbar
+        {194, -5132}, // Cascade+_bbar
+        {195, -5332} // Omega+_bbar
+      };
+
+      virtual bool generate(GenParticleCollection &);
+      void openFile(FILE *f);
+
+    private:
+
+      bool genEvent(std::map<std::pair<int,int>, GenParticleCollection> &particles_map);
+      float wrapvarBoxNo(const float var, const float low, const float high, int &boxno);
+      static float wrapvar(const float var, const float low, const float high);
+      std::vector<CLHEP::Hep3Vector> _targetBoxIntersections;
+      std::vector<CLHEP::Hep3Vector> _worldIntersections;
+      std::map<std::pair<int,int>, GenParticleCollection> _particles_map;
+
+      GlobalConstantsHandle<ParticleDataTable> pdt;
+
+      const float _GeV2MeV = CLHEP::GeV / CLHEP::MeV;
+      const float _ns2s = CLHEP::ns / CLHEP::s;
+      const float _cm2mm = CLHEP::cm / CLHEP::mm;
+
+      CLHEP::Hep3Vector _cosmicReferencePointInMu2e;
+      float _fluxConstant = 1.8e4;
+      float _tOffset = 0; ///< Time offset of sample, defaults to zero (no offset) [s]
+      bool _projectToTargetBox = false;
+      float _showerAreaExtension = 0; ///< Extend distribution of corsika particles in x,z by this much (e.g. 1000 will extend 10 m in -x, +x, -z, and +z) [mm]
+
+      std::string _refPointChoice;
+      bool _geomInfoObtained = false;
+      float _targetBoxXmin = 0;
+      float _targetBoxXmax = 0;
+      float _targetBoxYmin = 0;
+      float _targetBoxYmax = 0;
+      float _targetBoxZmin = 0;
+      float _targetBoxZmax = 0;
+
+      FILE *in = nullptr;
+      int _loops = 0; // number of loops, necessary to skip the garbage data between blocks
+      float _garbage;
+
+      unsigned int _primaries = 0;
+  };  // CosmicCORSIKA
+
+}
+
+#endif
\ No newline at end of file
diff --git a/Sources/src/CosmicCORSIKA.cc b/Sources/src/CosmicCORSIKA.cc
new file mode 100644
index 0000000000..ac68c0baea
--- /dev/null
+++ b/Sources/src/CosmicCORSIKA.cc
@@ -0,0 +1,259 @@
+////////////////////////////////////////////////////////////////////////
+/// \file  CosmicCORSIKA_module.cc
+/// \brief Generator for cosmic-ray secondaries based on pre-generated CORSIKA binary outputs.
+///
+/// \author Stefano Roberto Soleti roberto@lbl.gov
+////////////////////////////////////////////////////////////////////////
+
+#include "Sources/inc/CosmicCORSIKA.hh"
+
+using CLHEP::Hep3Vector;
+using CLHEP::HepLorentzVector;
+
+namespace mu2e {
+
+  CosmicCORSIKA::CosmicCORSIKA(const Config &conf)
+      : _fluxConstant(conf.fluxConstant()),
+        _tOffset(conf.tOffset()),
+        _projectToTargetBox(conf.projectToTargetBox()),
+        _showerAreaExtension(conf.showerAreaExtension()), // mm
+        _targetBoxXmin(conf.targetBoxXmin()), // mm
+        _targetBoxXmax(conf.targetBoxXmax()),  // mm
+        _targetBoxYmin(conf.targetBoxYmin()), // mm
+        _targetBoxYmax(conf.targetBoxYmax()),  // mm
+        _targetBoxZmin(conf.targetBoxZmin()), // mm
+        _targetBoxZmax(conf.targetBoxZmax())  // mm
+  {
+  }
+
+  const unsigned int CosmicCORSIKA::getNumShowers() {
+    return _primaries;
+  }
+
+  void CosmicCORSIKA::openFile(FILE *f) {
+    in = f;
+    fread(&_garbage, 4, 1, in);
+  }
+
+  const float CosmicCORSIKA::getLiveTime()
+  {
+    const float area = (_targetBoxXmax + 2 * _showerAreaExtension - _targetBoxXmin) * (_targetBoxZmax + 2 * _showerAreaExtension - _targetBoxZmin) * 1e-6; //m^2
+    const float eslope = -2.7;
+    const float lowE = 1.3;     // GeV
+    const float highE = 1e6; // GeV
+    const float EiToOneMinusGamma = pow(lowE, 1 + eslope);
+    const float EfToOneMinusGamma = pow(highE, 1 + eslope);
+    // http://pdg.lbl.gov/2018/reviews/rpp2018-rev-cosmic-rays.pdf eq. 29.2
+    return _primaries / (M_PI * area * _fluxConstant * (EfToOneMinusGamma - EiToOneMinusGamma) / (1. + eslope));
+  }
+
+  CosmicCORSIKA::~CosmicCORSIKA(){
+    std::cout << "Total number of primaries: " << getNumShowers() << std::endl;
+    std::cout << "Simulated live-time: " << getLiveTime() << std::endl;
+  }
+
+
+  float CosmicCORSIKA::wrapvarBoxNo(const float var, const float low, const float high, int &boxno)
+  {
+    //wrap variable so that it's always between low and high
+    boxno = int(floor(var / (high - low)));
+    return (var - (high - low) * floor(var / (high - low))) + low;
+  }
+
+  float CosmicCORSIKA::wrapvar( const float var, const float low, const float high){
+    //wrap variable so that it's always between low and high
+    return (var - (high - low) * floor(var/(high-low))) + low;
+  }
+
+
+ bool CosmicCORSIKA::genEvent(std::map<std::pair<int,int>, GenParticleCollection> &particles_map) {
+    float block[273]; // all data blocks have this size
+    // static TDatabasePDG *pdgt = TDatabasePDG::Instance();
+
+    bool running = true; // end run condition
+    bool validParticleSubBlock = true;
+    int particleDataSubBlockCounter = 0;
+
+    if (feof(in))
+      return false;
+
+    while (running)
+    {
+      std::vector<GenParticle> showerParticles;
+      fread(block, sizeof(block), 1, in); // blocks have 273 informations
+      char blockName[5];
+      memcpy(blockName, block, 4);
+      blockName[4] = '\0';
+      _loops++;
+
+      if (_loops % 21 == 0)
+      { // 2 _garbage data floats between every 21 blocks
+
+        fread(&_garbage, 4, 1, in);
+        fread(&_garbage, 4, 1, in);
+      }
+
+      // std::cout << blockName << " " << sizeof(block) << " " << (int)block[0] << std::endl;
+
+      //__________RUN HEADER
+
+      if (strcmp(blockName, "RUNH") == 0)
+      {
+        continue;
+      }
+
+      //__________EVENT HEADER
+      else if (strcmp(blockName, "EVTH") == 0)
+      {
+        _primaries++;
+        validParticleSubBlock = true;
+      }
+      //__________EVENT END
+
+      else if (strcmp(blockName, "EVTE") == 0)
+      {
+        if (particleDataSubBlockCounter > 0)
+        {
+          running = false;
+        }
+        particleDataSubBlockCounter = 0;
+      }
+      //__________END OF RUN
+      else if (strcmp(blockName, "RUNE") == 0)
+      {
+        _loops = 0;
+        return false;
+        // end run condition
+      }
+      else if (strcmp(blockName, "LONG")==0)
+      {
+        continue;
+      }
+      //__________PARTICLE DATA SUB-BLOCKS
+
+      else
+      {
+
+        if (validParticleSubBlock == true)
+        {
+          for (int l = 1; l < 40; l++)
+          { // each sub-block has up to 39 particles. If a sub-block is not fulfilled, trailing zeros are added
+
+            int k = 7 * (l - 1);
+            float id = block[k];
+            //cout << "	" << k << " " << block[k] << " " << block[k+1] << " " << block[k+2] << " " << block[k+3] << " " << block[k+4] << " " << block[k+5] << " " << block[k+6] << endl;	// for verification purposes only
+
+            if ((int)(id / 1000) != 0)
+            {
+              if (id < 75000)
+              {
+                // float dxyz[3] = {-block[k + 1], -block[k + 3], block[k + 2]};
+                // float xyz[3] = {-block[k + 4], FDHalfHeight, block[k + 5]};
+                // if (CheckTPCIntersection(xyz, dxyz)) {
+                if (block[k + 1] != 0)
+                {
+                  id = (int)(id / 1000);
+                  const int pdgId = corsikaToPdgId.at(id);
+                  const float P_x = block[k + 2] * _GeV2MeV;
+                  const float P_y = -block[k + 3] * _GeV2MeV;
+                  const float P_z = block[k + 1] * _GeV2MeV;
+
+                  int boxnox = 0, boxnoz = 0;
+
+                  const float x = wrapvarBoxNo(block[k + 5] * _cm2mm, _targetBoxXmin - _showerAreaExtension, _targetBoxXmax + _showerAreaExtension, boxnox);
+                  const float z = wrapvarBoxNo(-block[k + 4] * _cm2mm, _targetBoxZmin - _showerAreaExtension, _targetBoxZmax + _showerAreaExtension, boxnoz);
+                  std::pair xz(boxnox, boxnoz);
+                  const float m = pdt->particle(pdgId).ref().mass(); // to MeV
+
+                  const float energy = safeSqrt(P_x * P_x + P_y * P_y + P_z * P_z + m * m);
+
+                  const Hep3Vector position(x, 0, z);
+                  const HepLorentzVector mom4(P_x, P_y, P_z, energy);
+
+                  const float particleTime = block[k + 6];
+
+                  GenParticle part(static_cast<PDGCode::type>(pdgId),
+                                   GenId::cosmicCORSIKA, position, mom4,
+                                   particleTime);
+
+                  if (particles_map.count(xz) == 0) {
+                    GenParticleCollection parts;
+                    parts.push_back(part);
+                    particles_map.insert({xz, parts});
+                  } else {
+                    particles_map[xz].push_back(part);
+                  }
+
+                  particleDataSubBlockCounter++;
+                }
+              }
+            }
+            else
+            {
+              validParticleSubBlock = false;
+            }
+          }
+        }
+      }
+    }
+
+    return true;
+  }
+
+
+  bool CosmicCORSIKA::generate( GenParticleCollection& genParts)
+  {
+
+    // loop over particles in the truth object
+
+    bool passed = false;
+    while (!passed) {
+
+      if (_particles_map.size() == 0)
+      {
+        if (!genEvent(_particles_map))
+        {
+          return false;
+        }
+      }
+
+      GenParticleCollection particles = _particles_map.begin()->second;
+      GenParticleCollection crossingParticles;
+
+      float timeOffset = std::numeric_limits<float>::max();
+
+      for (unsigned int i = 0; i < particles.size(); i++) {
+        GenParticle particle = particles[i];
+
+
+        _targetBoxIntersections.clear();
+        VectorVolume particleTarget(particle.position(), particle.momentum().vect(),
+                                    _targetBoxXmin, _targetBoxXmax,
+                                    _targetBoxYmin, _targetBoxYmax,
+                                    _targetBoxZmin, _targetBoxZmax);
+
+        particleTarget.calIntersections(_targetBoxIntersections);
+
+        if (_targetBoxIntersections.size() > 0 || !_projectToTargetBox)
+          crossingParticles.push_back(particle);
+
+        if (particle.time() < timeOffset)
+          timeOffset = particle.time();
+      }
+
+      for (unsigned int i = 0; i < crossingParticles.size(); i++) {
+          GenParticle part = crossingParticles[i];
+          genParts.push_back(GenParticle(part.pdgId(), part.generatorId(), part.position(), part.momentum(), part.time()-timeOffset+_tOffset));
+      }
+      _particles_map.erase(_particles_map.begin()->first);
+
+      if (genParts.size() != 0)
+        passed = true;
+
+    }
+
+    return true;
+
+  }
+
+}// end namespace
diff --git a/Sources/src/FromCorsikaBinary_source.cc b/Sources/src/FromCorsikaBinary_source.cc
new file mode 100644
index 0000000000..debfce2532
--- /dev/null
+++ b/Sources/src/FromCorsikaBinary_source.cc
@@ -0,0 +1,211 @@
+// Transfer into framework CORSIKA binary files
+//
+// Original author: Stefano Roberto Soleti, 2019
+
+#include <iostream>
+#include <fstream>
+#include <boost/utility.hpp>
+#include <cassert>
+#include <set>
+#include <string>
+
+#include "CLHEP/Vector/LorentzVector.h"
+#include "CLHEP/Vector/ThreeVector.h"
+
+#include "art/Framework/IO/Sources/Source.h"
+#include "art/Framework/Core/InputSourceMacros.h"
+#include "art/Framework/IO/Sources/SourceHelper.h"
+#include "art/Framework/Principal/RunPrincipal.h"
+#include "art/Framework/Principal/SubRunPrincipal.h"
+#include "art/Framework/Principal/EventPrincipal.h"
+#include "art/Framework/IO/Sources/put_product_in_principal.h"
+#include "art/Utilities/Globals.h" // FIXME-KJK: should not be necessary to use this
+#include "canvas/Persistency/Provenance/Timestamp.h"
+#include "canvas/Persistency/Provenance/RunID.h"
+#include "canvas/Persistency/Provenance/SubRunID.h"
+#include "canvas/Persistency/Provenance/EventID.h"
+#include "canvas/Persistency/Provenance/BranchType.h"
+#include "canvas/Persistency/Provenance/ProductID.h"
+#include "canvas/Persistency/Provenance/canonicalProductName.h"
+#include "art/Framework/Core/ModuleMacros.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+#include "MCDataProducts/inc/GenParticle.hh"
+#include "MCDataProducts/inc/GenParticleCollection.hh"
+
+#include "Sources/inc/CosmicCORSIKA.hh"
+#include "SeedService/inc/SeedService.hh"
+
+using CLHEP::Hep3Vector;
+using CLHEP::HepLorentzVector;
+
+using namespace std;
+
+namespace mu2e {
+
+    //================================================================
+    class CorsikaBinaryDetail : private boost::noncopyable {
+      std::string myModuleLabel_;
+      art::SourceHelper const& pm_;
+      unsigned runNumber_; // from ParSet
+      art::SubRunID lastSubRunID_;
+      std::set<art::SubRunID> seenSRIDs_;
+
+      std::string currentFileName_;
+      FILE *currentFile_ = nullptr;
+      float garbage;
+
+      unsigned currentSubRunNumber_; // from file
+      // A helper function used to manage the principals.
+      // This is boilerplate that does not change if you change the data products.
+      void managePrincipals ( int runNumber,
+                              int subRunNumber,
+                              int eventNumber,
+                              art::RunPrincipal*&    outR,
+                              art::SubRunPrincipal*& outSR,
+                              art::EventPrincipal*&  outE);
+      unsigned getSubRunNumber(const std::string& filename) const;
+
+      unsigned currentEventNumber_;
+
+      art::ProductID particlesPID_;
+
+    public:
+      CorsikaBinaryDetail(const Parameters &conf,
+                          art::ProductRegistryHelper &,
+                          const art::SourceHelper &);
+
+      void readFile(std::string const& filename, art::FileBlock*& fb);
+
+      bool readNext(art::RunPrincipal* const& inR,
+                    art::SubRunPrincipal* const& inSR,
+                    art::RunPrincipal*& outR,
+                    art::SubRunPrincipal*& outSR,
+                    art::EventPrincipal*& outE);
+
+      void closeCurrentFile();
+
+      CosmicCORSIKA _corsikaGen;
+
+    };
+
+    //----------------------------------------------------------------
+    CorsikaBinaryDetail::CorsikaBinaryDetail(const Parameters& conf,
+                                             art::ProductRegistryHelper& rh,
+                                             const art::SourceHelper& pm)
+      : myModuleLabel_("FromCorsikaBinary")
+      , pm_(pm)
+      , runNumber_(conf().runNumber())
+      , currentSubRunNumber_(-1U)
+      , currentEventNumber_(0)
+      , _corsikaGen(conf())
+    {
+      if(!art::RunID(runNumber_).isValid()) {
+        throw cet::exception("BADCONFIG", " FromCorsikaBinary: ")
+          << " fhicl::ParameterSet specifies an invalid runNumber = "<<runNumber_<<"\n";
+      }
+
+      rh.reconstitutes<mu2e::GenParticleCollection,art::InEvent>(myModuleLabel_);
+
+      // auto get_ProductID = [](auto const& typeLabel, auto const& processName) {
+      //   if (!typeLabel.hasEmulatedModule()) {
+      //     throw cet::exception("BADCONFIG", " FromCorsikaBinary: ")
+      //     << " Must provided emulated module name for reconstituted product.\n";
+      //   }
+      //   auto const canonical_product_name = art::canonicalProductName(typeLabel.friendlyClassName(),
+      //                                                                 typeLabel.emulatedModule(),
+      //                                                                 typeLabel.productInstanceName(),
+      //                                                                 processName);
+      //   return art::ProductID{canonical_product_name};
+      // };
+      // particlesPID_ = get_ProductID(gpcTypeLabel, art::Globals::instance()->processName());
+
+    }
+
+    //----------------------------------------------------------------
+    void CorsikaBinaryDetail::readFile(const std::string& filename, art::FileBlock*& fb) {
+
+      currentFileName_ = filename;
+      currentSubRunNumber_ = getSubRunNumber(filename);
+      currentEventNumber_ = 0;
+
+      currentFile_ = fopen(filename.c_str(), "r");
+      _corsikaGen.openFile(currentFile_);
+      fb = new art::FileBlock(art::FileFormatVersion(1, "CorsikaBinaryInput"), currentFileName_);
+    }
+
+    //----------------------------------------------------------------
+    unsigned CorsikaBinaryDetail::getSubRunNumber(const std::string& filename) const {
+      const std::string::size_type islash = filename.find("DAT") ;
+      const std::string basename = (islash == std::string::npos) ? filename : filename.substr(islash + 3);
+      unsigned sr(-1);
+      std::istringstream is(basename);
+      if(!(is>>sr)) {
+        throw cet::exception("BADINPUTS")<<"Expect an unsigned integer at the beginning of input file name, got "<<basename<<"\n";
+      }
+      return sr;
+    }
+
+    //----------------------------------------------------------------
+    void CorsikaBinaryDetail::closeCurrentFile() {
+      currentFileName_ = "";
+      fclose(currentFile_);
+    }
+
+    //----------------------------------------------------------------
+    bool CorsikaBinaryDetail::readNext(art::RunPrincipal* const& inR,
+                                       art::SubRunPrincipal* const& inSR,
+                                       art::RunPrincipal*& outR,
+                                       art::SubRunPrincipal*& outSR,
+                                       art::EventPrincipal*& outE)
+    {
+      std::unique_ptr<GenParticleCollection> particles(new GenParticleCollection());
+      bool still_data = (_corsikaGen.generate(*particles));
+      if (!still_data) {
+        return false;
+      }
+
+      managePrincipals(runNumber_, currentSubRunNumber_, ++currentEventNumber_, outR, outSR, outE);
+      art::put_product_in_principal(std::move(particles), *outE, myModuleLabel_);
+
+      return true;
+
+    } // readNext()
+
+
+  // Each time that we encounter a new run, a new subRun or a new event, we need to make a new principal
+  // of the appropriate type.  This code does not need to change as the number and type of data products changes.
+  void CorsikaBinaryDetail::managePrincipals ( int runNumber,
+                                          int subRunNumber,
+                                          int eventNumber,
+                                          art::RunPrincipal*&    outR,
+                                          art::SubRunPrincipal*& outSR,
+                                          art::EventPrincipal*&  outE){
+
+    art::Timestamp ts;
+
+    art::SubRunID newID(runNumber, subRunNumber);
+    // if(newID.runID() != lastSubRunID_.runID()) {
+    //   // art takes ownership of the object pointed to by outR and will delete it at the appropriate time.
+    //   outR = pm_.makeRunPrincipal(runNumber, ts);
+    // }
+
+    if(newID != lastSubRunID_) {
+      // art takes ownership of the object pointed to by outSR and will delete it at the appropriate time.
+      outR = pm_.makeRunPrincipal(runNumber, ts);
+      outSR = pm_.makeSubRunPrincipal(runNumber,
+                                      subRunNumber,
+                                      ts);
+    }
+    lastSubRunID_ = newID;
+
+    // art takes ownership of the object pointed to by outE and will delete it at the appropriate time.
+    outE = pm_.makeEventPrincipal(runNumber, subRunNumber, eventNumber, ts, false);
+
+  } // managePrincipals()
+    //----------------------------------------------------------------
+
+} // namespace mu2e
+
+typedef art::Source<mu2e::CorsikaBinaryDetail> FromCorsikaBinary;
+DEFINE_ART_INPUT_SOURCE(FromCorsikaBinary);
diff --git a/Sources/src/SConscript b/Sources/src/SConscript
index 6c14cd9d8d..cb81566bab 100644
--- a/Sources/src/SConscript
+++ b/Sources/src/SConscript
@@ -9,13 +9,40 @@ Import('mu2e_helper')
 
 helper=mu2e_helper(env);
 
-mainlib = helper.make_mainlib ( [ 'mu2e_GlobalConstantsService',
-                                  'mu2e_MCDataProducts',
-                                  'cetlib',
-                                  'cetlib_except',
-                                  'HepPDT',
-                                  'CLHEP',
-                                  'boost_system'
+mainlib = helper.make_mainlib(['mu2e_GlobalConstantsService',
+                               'mu2e_Mu2eUtilities',
+                               'mu2e_ConditionsService',
+                               'mu2e_GeometryService',
+                               'mu2e_CalorimeterGeom',
+                               'mu2e_ExtinctionMonitorFNAL_Geometry',
+                               'mu2e_ProductionTargetGeom',
+                               'mu2e_StoppingTargetGeom',
+                               'mu2e_MCDataProducts',
+                               'mu2e_RecoDataProducts',
+                               'mu2e_Mu2eInterfaces',
+                               'mu2e_GlobalConstantsService',
+                               'mu2e_DataProducts',
+                               'mu2e_ConfigTools',
+                               'mu2e_GeneralUtilities',
+                               'art_Persistency_Provenance',
+                               'art_Persistency_Common',
+                               'art_Framework_Services_Registry',
+                               'art_Framework_Services_Optional_RandomNumberGenerator_service',
+                               'art_root_io_TFileService_service',
+                               'art_root_io_tfile_support',
+                               'art_Framework_Principal',
+                               'art_Framework_Core',
+                               'art_Utilities',
+                               'canvas',
+                               'fhiclcpp',
+                               'MF_MessageLogger',
+                               'cetlib',
+                               'cetlib_except',
+                               'CLHEP',
+                               'HepPDT',
+                               'HepPID',
+                               'boost_system',
+                               'gsl',
                                 ] )
 
 helper.make_plugins( [ mainlib,
@@ -37,6 +64,11 @@ helper.make_plugins( [ mainlib,
                        'Core',
                        'boost_filesystem',
                        'boost_system',
+                       'mu2e_GeometryService',
+                       'mu2e_CalorimeterGeom',
+                       'mu2e_ExtinctionMonitorFNAL_Geometry',
+                       'mu2e_ProductionTargetGeom',
+                       'mu2e_StoppingTargetGeom'
 ] )
 
 # This tells emacs to view this file in python mode.
