CosmicPionGen_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Module to generate only pions from cosmic rays
3 /// \author brebel@fnal.gov
4 /// \date February 2013
5 ////////////////////////////////////////////////////////////////////////
6 
7 extern "C" {
8 #include <getopt.h>
9 }
10 #include <cmath> // for fabs() etc
11 #include <memory>
12 
13 // ROOT includes
14 #include "TH1.h"
15 #include "TH2.h"
16 
17 // Framework includes
22 #include "fhiclcpp/ParameterSet.h"
28 
29 // NOvA includes
32 #include "Geometry/Geometry.h"
33 #include "NovaDAQConventions/DAQConventions.h"
35 #include "SummaryData/RunData.h"
36 #include "SummaryData/SubRunData.h"
39 
40 
41 namespace evgen {
42 
43  /// A module to check the results from the Monte Carlo generator
44  class CosmicPionGen : public art::EDProducer {
45  public:
46  explicit CosmicPionGen(fhicl::ParameterSet const& pset);
47  virtual ~CosmicPionGen();
48 
49  void produce(art::Event& evt);
50  void beginRun(art::Run& run);
51  void beginSubRun(art::SubRun& run);
52  void endSubRun(art::SubRun& run);
53 
54  private:
55 
56  evgb::CRYHelper* fCRYHelp; ///< CRY generator object
57  std::string fGeoVersion; ///< gdml file containing detector geometry
58  double fSpillLength; ///< Time in seconds to add to \ref fExposure every spill
59  double fExposure; ///< Livetime this subrun, in seconds
60  int fCycle; ///< MC production cycle
61  };
62 }
63 
64 
65 namespace evgen{
66 
67  //____________________________________________________________________________
69  : fSpillLength(pset.get< double >("SpillLength"))
70  , fExposure (0)
71  , fCycle (pset.get< int >("Cycle", 0))
72  {
73 
74  // Create a Art Random Number engine
75  int seed = pset.get< int >("Seed", evgb::GetRandomNumberSeed());
76  createEngine(seed);
77 
78  // get the random number generator service and make some CLHEP generators
80  CLHEP::HepRandomEngine& engine = rng->getEngine();
81 
82  // Create a new CRYHelper
83  fCRYHelp = new evgb::CRYHelper(pset, engine);
84 
85  // Producers
86  produces< std::vector<simb::MCTruth> >();
87  produces< sumdata::RunData, art::InRun >();
88  produces< sumdata::CosmicExposure, art::InSubRun >();
89  produces< sumdata::SubRunData, art::InSubRun >();
90  }
91 
92  //____________________________________________________________________________
94  {
95  if(fCRYHelp) delete fCRYHelp;
96  }
97 
98  //____________________________________________________________________________
100  {
101 
102  // grab the geometry object to see what geometry we are using
104 
105  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetId(),
106  geo->FileBaseName(),
107  geo->ExtractGDML()));
108 
109  run.put(std::move(runcol));
110 
111  return;
112  }
113 
114  //____________________________________________________________________________
116  {
117  fExposure = 0;
118  }
119 
120  //____________________________________________________________________________
122  {
123  std::unique_ptr<sumdata::CosmicExposure> ce(new sumdata::CosmicExposure);
124  ce->totlivetime = fExposure;
125  subrun.put(std::move(ce));
126 
127  // store the cycle information
129  std::unique_ptr< sumdata::SubRunData > sd(new sumdata::SubRunData(fCycle));
130  subrun.put(std::move(sd));
131 
132  }
133 
134  //____________________________________________________________________________
136  {
137  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
138 
140 
141  simb::MCTruth pitruth;
142  pitruth.SetOrigin(simb::kCosmicRay);
143 
144  while( pitruth.NParticles() < 1 ){
145  simb::MCTruth truth;
147  fCRYHelp->Sample(truth, geo->SurfaceY(), geo->DetLength(), 0);
148 
149  // look to see if there are any charged pions in this spill
150  // drop everything else and only write the pions to the event
151  for(int i = 0; i < truth.NParticles(); ++i){
153  if( abs(part.PdgCode()) == 211 || part.PdgCode() == 111 ) pitruth.Add(part);
154  }
155 
156  // add the spill length to the exposure even if
157  // we didn't have a pion so that we have the correct
158  // exposure in order to understand how long we would
159  // have to run to have a reasonable sample
161  }
162 
163  truthcol->push_back(pitruth);
164  evt.put(std::move(truthcol));
165 
166  return;
167  }
168 
169 }//namespace
170 
171 ////////////////////////////////////////////////////////////////////////
172 namespace evgen
173 {
175 }
Interface to the CRY cosmic-ray generator.
Definition: CRYHelper.h:26
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
unsigned int GetRandomNumberSeed()
Definition: evgenbase.h:22
base_engine_t & createEngine(seed_t seed)
double DetLength() const
void abs(TH1 *hist)
int NParticles() const
Definition: MCTruth.h:74
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
Definition: Run.h:31
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double sd(Eigen::VectorXd x)
base_engine_t & getEngine() const
TString part[npart]
Definition: Style.C:32
double Sample(simb::MCTruth &mctruth, double const &surfaceY, double const &detectorLength, double *w, double rantime=0)
Definition: CRYHelper.cxx:97
unsigned int seed
Definition: runWimpSim.h:102
Interface to the CRY cosmic ray generator.
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
void endSubRun(art::SubRun &run)
CosmicPionGen(fhicl::ParameterSet const &pset)
Definition: run.py:1
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
double fSpillLength
Time in seconds to add to fExposure every spill.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string fGeoVersion
gdml file containing detector geometry
void beginRun(art::Run &run)
double fExposure
Livetime this subrun, in seconds.
double SurfaceY() const
ProductID put(std::unique_ptr< PROD > &&)
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
Event generator information.
Definition: MCTruth.h:32
A module to check the results from the Monte Carlo generator.
Helper for AttenCurve.
Definition: Path.h:10
Module to generate only pions from cosmic rays.
void beginSubRun(art::SubRun &run)
Encapsulate the geometry of one entire detector (near, far, ndos)
Cosmic rays.
Definition: MCTruth.h:24
evgb::CRYHelper * fCRYHelp
CRY generator object.
int fCycle
MC production cycle.
void produce(art::Event &evt)
std::string FileBaseName() const