CosmicPionGen_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief  Module to generate only pions from cosmic rays
/// \author brebel@fnal.gov
/// \date   February 2013
////////////////////////////////////////////////////////////////////////

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

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

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

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


namespace evgen {

  /// A module to check the results from the Monte Carlo generator
  class CosmicPionGen : public art::EDProducer {
  public:
    explicit CosmicPionGen(fhicl::ParameterSet const& pset);
    virtual ~CosmicPionGen();                        

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

  private:

    evgb::CRYHelper* fCRYHelp;     ///< CRY generator object                    
    std::string      fGeoVersion;  ///< gdml file containing detector geometry
    double           fSpillLength; ///< Time in seconds to add to \ref fExposure every spill
    double           fExposure;    ///< Livetime this subrun, in seconds
    int              fCycle;       ///< MC production cycle
  };
}


namespace evgen{

  //____________________________________________________________________________
  CosmicPionGen::CosmicPionGen(fhicl::ParameterSet const& pset)
  : fSpillLength(pset.get< double >("SpillLength"))
  , fExposure   (0)
  , fCycle      (pset.get< int    >("Cycle", 0))
  {
      
    // Create a Art Random Number engine
    int seed = pset.get< int >("Seed", evgb::GetRandomNumberSeed());
    createEngine(seed);
    
    // get the random number generator service and make some CLHEP generators
    art::ServiceHandle<art::RandomNumberGenerator> rng;
    CLHEP::HepRandomEngine& engine = rng->getEngine();
    
    // Create a new CRYHelper
    fCRYHelp = new evgb::CRYHelper(pset, engine);
    
    // Producers
    produces< std::vector<simb::MCTruth> >();
    produces< sumdata::RunData,        art::InRun    >();
    produces< sumdata::CosmicExposure, art::InSubRun >();
    produces< sumdata::SubRunData,     art::InSubRun >();
  }
  
  //____________________________________________________________________________
  CosmicPionGen::~CosmicPionGen()
  {
    if(fCRYHelp) delete fCRYHelp;
  }

  //____________________________________________________________________________
  void CosmicPionGen::beginRun(art::Run& run)
  {
    
    // grab the geometry object to see what geometry we are using
    art::ServiceHandle<geo::Geometry> geo;
    
    std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetId(),
                                                                  geo->FileBaseName(),
                                                                  geo->ExtractGDML()));
    
    run.put(std::move(runcol));

    return;
  }

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

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

  }

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

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

    simb::MCTruth pitruth;
    pitruth.SetOrigin(simb::kCosmicRay);

    while( pitruth.NParticles() < 1 ){
      simb::MCTruth  truth;
      truth.SetOrigin(simb::kCosmicRay);
      fCRYHelp->Sample(truth, geo->SurfaceY(), geo->DetLength(), 0);

      // look to see if there are any charged pions in this spill
      // drop everything else and only write the pions to the event
      for(int i = 0; i < truth.NParticles(); ++i){
        simb::MCParticle part(truth.GetParticle(i));
        if( abs(part.PdgCode()) == 211 || part.PdgCode() == 111 ) pitruth.Add(part);
      }

      // add the spill length to the exposure even if 
      // we didn't have a pion so that we have the correct
      // exposure in order to understand how long we would
      // have to run to have a reasonable sample
      fExposure += fSpillLength;
    }

    truthcol->push_back(pitruth);
    evt.put(std::move(truthcol));

    return;
  }

}//namespace

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