LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
evgen::CORSIKAGen Class Reference

LArSoft interface to CORSIKA event generator. More...

Inheritance diagram for evgen::CORSIKAGen:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 CORSIKAGen (fhicl::ParameterSet const &pset)
 
virtual ~CORSIKAGen ()
 
void produce (art::Event &evt)
 
void beginRun (art::Run &run)
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
void fillProductDescriptions ()
 
void registerProducts (ProductDescriptions &productsToRegister)
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
std::unique_ptr< Worker > makeWorker (WorkerParams const &wp)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Protected Member Functions

ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Private Member Functions

void openDBs (std::string const &module_label)
 
void populateNShowers ()
 
void populateTOffset ()
 
void GetSample (simb::MCTruth &)
 
double wrapvar (const double var, const double low, const double high)
 
double wrapvarBoxNo (const double var, const double low, const double high, int &boxno)
 
void ProjectToBoxEdge (const double xyz[], const double dxyz[], const double xlo, const double xhi, const double ylo, const double yhi, const double zlo, const double zhi, double xyzout[])
 Propagates a point back to the surface of a box. More...
 

Private Attributes

int fShowerInputs = 0
 Number of shower inputs to process from. More...
 
std::vector< double > fNShowersPerEvent
 Number of showers to put in each event of duration fSampleTime; one per showerinput. More...
 
std::vector< int > fMaxShowers
 
double fShowerBounds [6]
 Boundaries of area over which showers are to be distributed (x(min), x(max), unused, y, z(min), z(max) ) More...
 
double fToffset_corsika
 Timing offset to account for propagation time through atmosphere, populated from db. More...
 
ifdh_ns::ifdh * fIFDH = 0
 (optional) flux file handling More...
 
double fProjectToHeight = 0.
 Height to which particles will be projected [cm]. More...
 
std::vector< std::string > fShowerInputFiles
 Set of CORSIKA shower data files to use. More...
 
std::string fShowerCopyType
 
std::vector< double > fShowerFluxConstants
 Set of flux constants to be associated with each shower data file. More...
 
double fSampleTime = 0.
 Duration of sample [s]. More...
 
double fToffset = 0.
 Time offset of sample, defaults to zero (no offset) [s]. More...
 
std::vector< double > fBuffBox
 Buffer box extensions to cryostat in each direction (6 of them: x_lo,x_hi,y_lo,y_hi,z_lo,z_hi) [cm]. More...
 
double fShowerAreaExtension
 Extend distribution of corsika particles in x,z by this much (e.g. 1000 will extend 10 m in -x, +x, -z, and +z) [cm]. More...
 
sqlite3 * fdb [5]
 Pointers to sqlite3 database object, max of 5. More...
 
double fRandomXZShift
 Each shower will be shifted by a random amount in xz so that showers won't repeatedly sample the same space [cm]. More...
 
CLHEP::HepRandomEngine & fGenEngine
 
CLHEP::HepRandomEngine & fPoisEngine
 

Detailed Description

LArSoft interface to CORSIKA event generator.

Note
A presentation on this module by the original author is archived at: https://indico.fnal.gov/event/10893/contribution/3/material/slides

In CORSIKA jargon, a "shower" is the cascade of particles resulting from a primary cosmic ray interaction. This module creates a single simb::MCTruth object (stored as data product into a std::vector<simb::MCTruth> with a single entry) containing all the particles from cosmic ray showers crossing a surface above the detector.

The generation procedure consists of selecting showers from a database of pregenerated events, and then to adapt them to the parameters requested in the module configuration. Pregenerated showers are "observed" at a altitude set in CORSIKA configuration.

Databases need to be stored as files in SQLite3 format. Multiple file sources can be specified (ShowerInputFiles configuration parameter). From each source, one database file is selected and copied locally via IFDH. From each source, showers are extracted proportionally to the relative flux specified in the configuration (specified in ShowerFluxConstants, see normalization below). The actual number of showers per event and per source is extracted according to a Poisson distribution around the predicted average number of primary cosmic rays for that source.

Flux normalization

CORSIKA generates showers from each specific cosmic ray type $ A $ (e.g. iron, proton, etc.) according to a power law distribution $ \Phi_{A}(E) \propto E^{-\gamma_{A}} $ of the primary particle energy $ E $ [GeV]. When sampling pregenerated events, we bypass the normalization imposed by CORSIKA and gain complete control on it.

Within CORSIKAGen, for each source (usually each on a different primary cosmic ray type, e.g. iron, proton, etc.), the average number of generated showers is $ n_{A} = \pi S T k_{A} \int E^{-\gamma_{A}} dE $ with $ S $ the area of the surface the flux passes across, $ T $ the exposure time, the integral defined over the full energy range of the pregenerated showers in the source, and $ k_{A} $ a factor specified in the configuration (ShowerFluxConstants parameters). This is the flux of primary cosmic rays, not of the observed particles from their showers. Note that it depends on an area and a time interval, but it is uniform with respect to translations and constant in time.

As explained below, we consider only the secondary particles that cross an "observation" surface. After cosmic ray primary particles cross the flux surface ( $ S_{\Phi} $ above), they develop into showers of particles that spread across large areas. Limiting ourself to the observation of particles on a small surface $ S_{o} $ has two effects. We lose the part of the showers that misses that surface $ S_{o} $. Also, considering a span of time with multiple showers, we miss particles from other hypothetical showers whose primaries crossed outside the generation surface $ S_{\Phi} $ whose shower development would leak into $ S_{o} $: we did not simulate those showers at all! In terms of total flux of observed particles, under the assumption that the flux is uniform in space, choosing $ S_{\Phi} $ the same size as $ S_{o} $ makes the two effects get the same size: just as many particles from primaries from $ S_{\Phi} $ leak out of $ S_{o} $, as many particles from primaries from outside $ S_{\Phi} $ sneak in $ S_{o} $. In that case, counting all the particles from the primaries crossing a surface $ S_{\Phi} $ of area S regardless of where they land yields the right amount of cosmic ray secondary particles across an observation surface $ S_{o} $ also of area S. Practically, the particles landing outside $ S_{o} $ need to be recovered to preserve the correct normalization, which is described in the next section.

Surface coverage, position and timing

The surface we detect the particles through (let's call it $ S_{d} $) is defined by the smallest rectangle including all cryostats in the detector, and located at the height of the ceiling of the tallest cryostat. This surface can be increased by specifying a positive value for ShowerAreaExtension configuration parameter, in which case each side of the rectangle will be extended by that amount.

Showers are extracted one by one from the pregenerated samples and treated independently. Ideally, the detection surface $ S_{d} $ would be at the same exact altitude as the observation surface set in CORSIKA (called $ S_{o} $ above). In practice, we go the other way around, with the assumption that the shower observed at $ S_{d} $ would be very close to the one we actually generated at $ S_{o} $, and teleport the generated particles on $ S_{d} $. Since the cryostats may be just meters from the earth surface $ S_{o} $ lies on, this is an acceptable approximation.

All the particles of one shower are compelled within surface $ S_{d} $ as a first step. As explained when describing the "normalization", we need to keep all the shower particles, one way or the other. So, particles of the shower that fell out of $ S_{d} $ are repackaged into other showers and translated back in. This is equivalent to assume the primary cosmic rays originating such shower would happen outside our generation volume ( $ S_{\Phi} $), and their shower would then spill into $ S_{d} $. Since these repackaged showers are in principle independent and uncorrelated, they are assigned a random time different than the main shower, leveraging the assumption of constantness of the flux.

As for the azimuth, this module uses an approximation by setting north direction to match the z axis of LArSoft geometry (typically assumed to be the direction of the beam particle).

The particles so manipulated are then back-propagated from the observation surface to an absolute height defined by ProjectToHeight (although for particular combination of position and direction, the particles might be propagated back to the edge of the world, or even outside it).

As a final filter, only the particles whose straight projections cross any of the cryostats (with some buffer volume around, defined by BufferBox) are stored, while the other ones are discarded. Note that the actual interactions that particles from the generation surface undergo may deviate them enough to miss the cryostats anyway, and conversely particles that have been filtered out because shooting off the cryostats might have been subsequently deviated to actually cross them. This effect is not corrected for at this time.

The time of the showers is uniformly distributed within the configured time interval, defined by SampleTime starting from TimeOffset.

Configuration parameters

  • ShowerInputFiles (list of paths; mandatory): a list of file paths to pregenerated CORSIKA shower files. Each entry can be a single file or use wildcards (*) to specify a set of files to choose among. Paths and wildcards are processed by IFDH.
  • ShowerFluxConstants (list of real numbers; mandatory): for each entry $ A $ in ShowerInputFiles, specify the normalization factor $ K_{A} $ of their distribution [ $ m^{-2}s^{-1} $]
  • ProjectToHeight (real, default: 0): the generated particles will appear to come from this height [cm]
  • TimeOffset (real; default: 0): start time of the exposure window [s], relative to the simulation time start
  • SampleTime (real; mandatory): duration of the simulated exposure to cosmic rays [s]
  • ShowerAreaExtension (real; default: 0): extend the size of the observation surface of shower particles (S) by this much [cm]; e.g. 1000 will extend 10 m on each side
  • RandomXZShift (real; default: 0): the original position of each shower is randomly shifted within a square with this length as side [cm]
  • BufferBox (list of six lengths, all 0 by default): extension to the volume of each cryostat for the purpose of filtering out the particles which do not cross the detector; each cryostat volume is independently extended by the same amount, specified here as shifts to lower x, higher x, lower y, higher y, lower z and higher z, in that order cm
  • SeedGenerator (integer): force random number generator for event generation to the specified value
  • SeedPoisson (integer): force random number generator for number of showers to the specified value
  • Seed: alias for SeedGenerator

Random engines

Currently two random engines are used:

  • a generator engine (driven by SeedGenerator), of general use
  • a "Poisson" engine (driven by SeedPoisson), only used to determine the number of showers to be selected on each event

Definition at line 221 of file CORSIKAGen_module.cc.

Member Typedef Documentation

Definition at line 17 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 26 of file Producer.h.

Constructor & Destructor Documentation

evgen::CORSIKAGen::CORSIKAGen ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 302 of file CORSIKAGen_module.cc.

References art::detail::EngineCreator::createEngine(), fBuffBox, fGenEngine, fPoisEngine, fProjectToHeight, fRandomXZShift, fSampleTime, fShowerAreaExtension, fShowerCopyType, fShowerFluxConstants, fShowerInputFiles, fShowerInputs, fToffset, openDBs(), populateNShowers(), and populateTOffset().

303  : EDProducer{p}
304  , fProjectToHeight(p.get<double>("ProjectToHeight", 0.))
305  , fShowerInputFiles(p.get<std::vector<std::string>>("ShowerInputFiles"))
306  , fShowerCopyType(p.get<std::string>("ShowerCopyType", "IFDH"))
307  , fShowerFluxConstants(p.get<std::vector<double>>("ShowerFluxConstants"))
308  , fSampleTime(p.get<double>("SampleTime", 0.))
309  , fToffset(p.get<double>("TimeOffset", 0.))
310  , fBuffBox(p.get<std::vector<double>>("BufferBox", {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}))
311  , fShowerAreaExtension(p.get<double>("ShowerAreaExtension", 0.))
312  , fRandomXZShift(p.get<double>("RandomXZShift", 0.))
313  , fGenEngine(art::ServiceHandle<rndm::NuRandomService>()->registerAndSeedEngine(
314  createEngine(0, "HepJamesRandom", "gen"),
315  "HepJamesRandom",
316  "gen",
317  p,
318  {"Seed", "SeedGenerator"}))
319  , fPoisEngine(art::ServiceHandle<rndm::NuRandomService>()->registerAndSeedEngine(
320  createEngine(0, "HepJamesRandom", "pois"),
321  "HepJamesRandom",
322  "pois",
323  p,
324  "SeedPoisson"))
325  {
326  if (fShowerInputFiles.size() != fShowerFluxConstants.size() || fShowerInputFiles.size() == 0 ||
327  fShowerFluxConstants.size() == 0)
328  throw cet::exception("CORSIKAGen")
329  << "ShowerInputFiles and ShowerFluxConstants have different or invalid sizes!"
330  << "\n";
332 
333  if (fSampleTime == 0.) throw cet::exception("CORSIKAGen") << "SampleTime not set!";
334 
335  if (fProjectToHeight == 0.) mf::LogInfo("CORSIKAGen") << "Using 0. for fProjectToHeight!";
336  // create a default random engine; obtain the random seed from NuRandomService,
337  // unless overridden in configuration with key "Seed"
338 
339  this->openDBs(p.get<std::string>("module_label"));
340  this->populateNShowers();
341  this->populateTOffset();
342 
343  produces<std::vector<simb::MCTruth>>();
344  produces<sumdata::RunData, art::InRun>();
345  }
base_engine_t & createEngine(seed_t seed)
double fToffset
Time offset of sample, defaults to zero (no offset) [s].
double fSampleTime
Duration of sample [s].
double fProjectToHeight
Height to which particles will be projected [cm].
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< std::string > fShowerInputFiles
Set of CORSIKA shower data files to use.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
int fShowerInputs
Number of shower inputs to process from.
CLHEP::HepRandomEngine & fPoisEngine
std::string fShowerCopyType
std::vector< double > fShowerFluxConstants
Set of flux constants to be associated with each shower data file.
double fRandomXZShift
Each shower will be shifted by a random amount in xz so that showers won&#39;t repeatedly sample the same...
CLHEP::HepRandomEngine & fGenEngine
void openDBs(std::string const &module_label)
double fShowerAreaExtension
Extend distribution of corsika particles in x,z by this much (e.g. 1000 will extend 10 m in -x...
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< double > fBuffBox
Buffer box extensions to cryostat in each direction (6 of them: x_lo,x_hi,y_lo,y_hi,z_lo,z_hi) [cm].
evgen::CORSIKAGen::~CORSIKAGen ( )
virtual

Definition at line 347 of file CORSIKAGen_module.cc.

References fdb, fIFDH, and fShowerInputs.

348  {
349  for (int i = 0; i < fShowerInputs; i++) {
350  sqlite3_close(fdb[i]);
351  }
352  //cleanup temp files
353  fIFDH->cleanup();
354  }
ifdh_ns::ifdh * fIFDH
(optional) flux file handling
int fShowerInputs
Number of shower inputs to process from.
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.

Member Function Documentation

void evgen::CORSIKAGen::beginRun ( art::Run run)
virtual

Reimplemented from art::EDProducer.

Definition at line 801 of file CORSIKAGen_module.cc.

References geo::GeometryCore::DetectorName(), art::fullRun(), and art::Run::put().

802  {
804  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
805  }
constexpr auto fullRun()
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Run.h:121
std::string const & DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:203
Namespace collecting geometry-related classes utilities.
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::consumes ( InputTag const &  tag)
protectedinherited

Definition at line 61 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumes().

62  {
63  return collector_.consumes<T, BT>(tag);
64  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
ProductToken< T > consumes(InputTag const &)
ConsumesCollector & art::ModuleBase::consumesCollector ( )
protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

58  {
59  return collector_;
60  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename T , BranchType BT>
void art::ModuleBase::consumesMany ( )
protectedinherited

Definition at line 75 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesMany().

76  {
77  collector_.consumesMany<T, BT>();
78  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::consumesView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::consumesView ( InputTag const &  tag)
inherited

Definition at line 68 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesView().

69  {
70  return collector_.consumesView<T, BT>(tag);
71  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
ViewToken< Element > consumesView(InputTag const &)
void art::detail::Producer::doBeginJob ( SharedResources const &  resources)
inherited

Definition at line 22 of file Producer.cc.

References art::detail::Producer::beginJobWithFrame(), and art::detail::Producer::setupQueues().

23  {
24  setupQueues(resources);
25  ProcessingFrame const frame{ScheduleID{}};
26  beginJobWithFrame(frame);
27  }
virtual void setupQueues(SharedResources const &)=0
virtual void beginJobWithFrame(ProcessingFrame const &)=0
bool art::detail::Producer::doBeginRun ( RunPrincipal rp,
ModuleContext const &  mc 
)
inherited

Definition at line 65 of file Producer.cc.

References art::detail::Producer::beginRunWithFrame(), art::RangeSet::forRun(), art::RunPrincipal::makeRun(), r, art::RunPrincipal::runID(), and art::ModuleContext::scheduleID().

66  {
67  auto r = rp.makeRun(mc, RangeSet::forRun(rp.runID()));
68  ProcessingFrame const frame{mc.scheduleID()};
69  beginRunWithFrame(r, frame);
70  r.commitProducts();
71  return true;
72  }
TRandom r
Definition: spectrum.C:23
virtual void beginRunWithFrame(Run &, ProcessingFrame const &)=0
static RangeSet forRun(RunID)
Definition: RangeSet.cc:51
bool art::detail::Producer::doBeginSubRun ( SubRunPrincipal srp,
ModuleContext const &  mc 
)
inherited

Definition at line 85 of file Producer.cc.

References art::detail::Producer::beginSubRunWithFrame(), art::RangeSet::forSubRun(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::SubRunPrincipal::subRunID().

86  {
87  auto sr = srp.makeSubRun(mc, RangeSet::forSubRun(srp.subRunID()));
88  ProcessingFrame const frame{mc.scheduleID()};
89  beginSubRunWithFrame(sr, frame);
90  sr.commitProducts();
91  return true;
92  }
virtual void beginSubRunWithFrame(SubRun &, ProcessingFrame const &)=0
static RangeSet forSubRun(SubRunID)
Definition: RangeSet.cc:57
void art::detail::Producer::doEndJob ( )
inherited

Definition at line 30 of file Producer.cc.

References art::detail::Producer::endJobWithFrame().

31  {
32  ProcessingFrame const frame{ScheduleID{}};
33  endJobWithFrame(frame);
34  }
virtual void endJobWithFrame(ProcessingFrame const &)=0
bool art::detail::Producer::doEndRun ( RunPrincipal rp,
ModuleContext const &  mc 
)
inherited

Definition at line 75 of file Producer.cc.

References art::detail::Producer::endRunWithFrame(), art::RunPrincipal::makeRun(), r, art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

76  {
77  auto r = rp.makeRun(mc, rp.seenRanges());
78  ProcessingFrame const frame{mc.scheduleID()};
79  endRunWithFrame(r, frame);
80  r.commitProducts();
81  return true;
82  }
TRandom r
Definition: spectrum.C:23
virtual void endRunWithFrame(Run &, ProcessingFrame const &)=0
bool art::detail::Producer::doEndSubRun ( SubRunPrincipal srp,
ModuleContext const &  mc 
)
inherited

Definition at line 95 of file Producer.cc.

References art::detail::Producer::endSubRunWithFrame(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

96  {
97  auto sr = srp.makeSubRun(mc, srp.seenRanges());
98  ProcessingFrame const frame{mc.scheduleID()};
99  endSubRunWithFrame(sr, frame);
100  sr.commitProducts();
101  return true;
102  }
virtual void endSubRunWithFrame(SubRun &, ProcessingFrame const &)=0
bool art::detail::Producer::doEvent ( EventPrincipal ep,
ModuleContext const &  mc,
std::atomic< std::size_t > &  counts_run,
std::atomic< std::size_t > &  counts_passed,
std::atomic< std::size_t > &  counts_failed 
)
inherited

Definition at line 105 of file Producer.cc.

References art::detail::Producer::checkPutProducts_, e, art::EventPrincipal::makeEvent(), art::detail::Producer::produceWithFrame(), and art::ModuleContext::scheduleID().

110  {
111  auto e = ep.makeEvent(mc);
112  ++counts_run;
113  ProcessingFrame const frame{mc.scheduleID()};
114  produceWithFrame(e, frame);
115  e.commitProducts(checkPutProducts_, &expectedProducts<InEvent>());
116  ++counts_passed;
117  return true;
118  }
bool const checkPutProducts_
Definition: Producer.h:70
Float_t e
Definition: plot.C:35
virtual void produceWithFrame(Event &, ProcessingFrame const &)=0
void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const &  fb)
inherited

Definition at line 44 of file Producer.cc.

References art::detail::Producer::respondToCloseInputFileWithFrame().

45  {
46  ProcessingFrame const frame{ScheduleID{}};
48  }
virtual void respondToCloseInputFileWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const &  fb)
inherited

Definition at line 58 of file Producer.cc.

References art::detail::Producer::respondToCloseOutputFilesWithFrame().

59  {
60  ProcessingFrame const frame{ScheduleID{}};
62  }
virtual void respondToCloseOutputFilesWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const &  fb)
inherited

Definition at line 37 of file Producer.cc.

References art::detail::Producer::respondToOpenInputFileWithFrame().

38  {
39  ProcessingFrame const frame{ScheduleID{}};
41  }
virtual void respondToOpenInputFileWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const &  fb)
inherited

Definition at line 51 of file Producer.cc.

References art::detail::Producer::respondToOpenOutputFilesWithFrame().

52  {
53  ProcessingFrame const frame{ScheduleID{}};
55  }
virtual void respondToOpenOutputFilesWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::Modifier::fillProductDescriptions ( )
inherited

Definition at line 10 of file Modifier.cc.

References art::ProductRegistryHelper::fillDescriptions(), and art::ModuleBase::moduleDescription().

11  {
13  }
void fillDescriptions(ModuleDescription const &md)
ModuleDescription const & moduleDescription() const
Definition: ModuleBase.cc:13
std::array< std::vector< ProductInfo >, NumBranchTypes > const & art::ModuleBase::getConsumables ( ) const
inherited

Definition at line 43 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::getConsumables().

44  {
45  return collector_.getConsumables();
46  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables() const
void evgen::CORSIKAGen::GetSample ( simb::MCTruth mctruth)
private

Definition at line 625 of file CORSIKAGen_module.cc.

References simb::MCTruth::Add(), simb::MCParticle::AddTrajectoryPoint(), fdb, fGenEngine, fMaxShowers, fNShowersPerEvent, fPoisEngine, fProjectToHeight, fRandomXZShift, fSampleTime, fShowerBounds, fShowerInputs, fToffset, fToffset_corsika, MF_LOG_DEBUG, ProjectToBoxEdge(), geo::GeometryCore::WorldBox(), wrapvar(), wrapvarBoxNo(), x, x1, x2, y1, y2, and z.

Referenced by produce().

626  {
627  //for each input, randomly pull fNShowersPerEvent[i] showers from the Particles table
628  //and randomly place them in time (between -fSampleTime/2 and fSampleTime/2)
629  //wrap their positions based on the size of the area under consideration
630  //based on http://nusoft.fnal.gov/larsoft/doxsvn/html/CRYHelper_8cxx_source.html (Sample)
631 
632  //query from sqlite db with select * from particles where shower in (select id from showers ORDER BY substr(id*0.51123124141,length(id)+2) limit 100000) ORDER BY substr(shower*0.51123124141,length(shower)+2);
633  //where 0.51123124141 is a random seed to allow randomly selecting rows and should be randomly generated for each query
634  //the inner order by is to select randomly from the possible shower id's
635  //the outer order by is to make sure the shower numbers are ordered randomly (without this, the showers always come out ordered by shower number
636  //and 100000 is the number of showers to be selected at random and needs to be less than the number of showers in the showers table
637 
638  //TDatabasePDG is for looking up particle masses
639  static TDatabasePDG* pdgt = TDatabasePDG::Instance();
640 
641  //db variables
642  sqlite3_stmt* statement;
643  const TString kStatement(
644  "select shower,pdg,px,py,pz,x,z,t,e from particles where shower in (select id from showers "
645  "ORDER BY substr(id*%f,length(id)+2) limit %d) ORDER BY substr(shower*%f,length(shower)+2)");
646 
647  CLHEP::RandFlat flat(fGenEngine);
648  CLHEP::RandPoissonQ randpois(fPoisEngine);
649 
650  // get geometry and figure where to project particles to, based on CRYHelper
652  double x1, x2;
653  double y1, y2;
654  double z1, z2;
655  geom->WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);
656 
657  // make the world box slightly smaller so that the projection to
658  // the edge avoids possible rounding errors later on with Geant4
659  double fBoxDelta = 1.e-5;
660  x1 += fBoxDelta;
661  x2 -= fBoxDelta;
662  y1 += fBoxDelta;
663  y2 = fProjectToHeight;
664  z1 += fBoxDelta;
665  z2 -= fBoxDelta;
666 
667  //populate mctruth
668  int ntotalCtr = 0; //count number of particles added to mctruth
669  int lastShower =
670  0; //keep track of last shower id so that t can be randomized on every new shower
671  int nShowerCntr = 0; //keep track of how many showers are left to be added to mctruth
672  int nShowerQry = 0; //number of showers to query from db
673  int shower, pdg;
674  double px, py, pz, x, z, tParticleTime, etot,
675  showerTime = 0., showerTimex = 0., showerTimez = 0., showerXOffset = 0., showerZOffset = 0.,
676  t;
677  for (int i = 0; i < fShowerInputs; i++) {
678  nShowerCntr = randpois.fire(fNShowersPerEvent[i]);
679  mf::LogInfo("CORSIKAGEN") << " Shower input " << i << " with mean " << fNShowersPerEvent[i]
680  << " generating " << nShowerCntr;
681 
682  while (nShowerCntr > 0) {
683  //how many showers should we query?
684  if (nShowerCntr > fMaxShowers[i]) {
685  nShowerQry = fMaxShowers[i]; //take the group size
686  }
687  else {
688  nShowerQry = nShowerCntr; //take the rest that are needed
689  }
690  //build and do query to get nshowers
691  double thisrnd = flat(); //need a new random number for each query
692  TString kthisStatement = TString::Format(kStatement.Data(), thisrnd, nShowerQry, thisrnd);
693  MF_LOG_DEBUG("CORSIKAGen") << "Executing: " << kthisStatement;
694  if (sqlite3_prepare(fdb[i], kthisStatement.Data(), -1, &statement, 0) == SQLITE_OK) {
695  int res = 0;
696  //loop over database rows, pushing particles into mctruth object
697  while (1) {
698  res = sqlite3_step(statement);
699  if (res == SQLITE_ROW) {
700  /*
701  * Memo columns:
702  * [0] shower
703  * [1] particle ID (PDG)
704  * [2] momentum: x component [GeV/c]
705  * [3] momentum: y component [GeV/c]
706  * [4] momentum: z component [GeV/c]
707  * [5] position: x component [cm]
708  * [6] position: z component [cm]
709  * [7] time [ns]
710  * [8] energy [GeV]
711  */
712  shower = sqlite3_column_int(statement, 0);
713  if (shower != lastShower) {
714  //each new shower gets its own random time and position offsets
715  showerTime = 1e9 * (flat() * fSampleTime); //converting from s to ns
716  showerTimex = 1e9 * (flat() * fSampleTime); //converting from s to ns
717  showerTimez = 1e9 * (flat() * fSampleTime); //converting from s to ns
718  //and a random offset in both z and x controlled by the fRandomXZShift parameter
719  showerXOffset = flat() * fRandomXZShift - (fRandomXZShift / 2);
720  showerZOffset = flat() * fRandomXZShift - (fRandomXZShift / 2);
721  }
722  pdg = sqlite3_column_int(statement, 1);
723  //get mass for this particle
724  double m = 0.; // in GeV
725  TParticlePDG* pdgp = pdgt->GetParticle(pdg);
726  if (pdgp) m = pdgp->Mass();
727 
728  //Note: position/momentum in db have north=-x and west=+z, rotate so that +z is north and +x is west
729  //get momentum components
730  px = sqlite3_column_double(statement, 4); //uboone x=Particlez
731  py = sqlite3_column_double(statement, 3);
732  pz = -sqlite3_column_double(statement, 2); //uboone z=-Particlex
733  etot = sqlite3_column_double(statement, 8);
734 
735  //get/calculate position components
736  int boxnoX = 0, boxnoZ = 0;
737  x = wrapvarBoxNo(sqlite3_column_double(statement, 6) + showerXOffset,
738  fShowerBounds[0],
739  fShowerBounds[1],
740  boxnoX);
741  z = wrapvarBoxNo(-sqlite3_column_double(statement, 5) + showerZOffset,
742  fShowerBounds[4],
743  fShowerBounds[5],
744  boxnoZ);
745  tParticleTime = sqlite3_column_double(
746  statement, 7); //time offset, includes propagation time from top of atmosphere
747  //actual particle time is particle surface arrival time
748  //+ shower start time
749  //+ global offset (fcl parameter, in s)
750  //- propagation time through atmosphere
751  //+ boxNo{X,Z} time offset to make grid boxes have different shower times
752  t = tParticleTime + showerTime + (1e9 * fToffset) - fToffset_corsika +
753  showerTimex * boxnoX + showerTimez * boxnoZ;
754  //wrap surface arrival so that it's in the desired time window
755  t = wrapvar(t, (1e9 * fToffset), 1e9 * (fToffset + fSampleTime));
756 
757  simb::MCParticle p(ntotalCtr, pdg, "primary", -200, m, 1);
758 
759  //project back to wordvol/fProjectToHeight
760  /*
761  * This back propagation goes from a point on the upper surface of
762  * the cryostat back to the edge of the world, except that that
763  * world is cut short by `fProjectToHeight` (`y2`) ceiling.
764  * The projection will most often lie on that ceiling, but it may
765  * end up instead on one of the side edges of the world, or even
766  * outside it.
767  */
768  double xyzo[3];
769  double x0[3] = {x, fShowerBounds[3], z};
770  double dx[3] = {px, py, pz};
771  this->ProjectToBoxEdge(x0, dx, x1, x2, y1, y2, z1, z2, xyzo);
772 
773  TLorentzVector pos(
774  xyzo[0], xyzo[1], xyzo[2], t); // time needs to be in ns to match GENIE, etc
775  TLorentzVector mom(px, py, pz, etot);
776  p.AddTrajectoryPoint(pos, mom);
777  mctruth.Add(p);
778  ntotalCtr++;
779  lastShower = shower;
780  }
781  else if (res == SQLITE_DONE) {
782  break;
783  }
784  else {
785  throw cet::exception("CORSIKAGen")
786  << "Unexpected sqlite3_step return value: (" << res << "); "
787  << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
788  }
789  }
790  }
791  else {
792  throw cet::exception("CORSIKAGen")
793  << "Error preparing statement: (" << kthisStatement << "); "
794  << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
795  }
796  nShowerCntr = nShowerCntr - nShowerQry;
797  }
798  }
799  }
Float_t x
Definition: compare.C:6
double fToffset
Time offset of sample, defaults to zero (no offset) [s].
double wrapvar(const double var, const double low, const double high)
double fSampleTime
Duration of sample [s].
double fProjectToHeight
Height to which particles will be projected [cm].
Float_t y1[n_points_granero]
Definition: compare.C:5
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void WorldBox(double *xlo, double *xhi, double *ylo, double *yhi, double *zlo, double *zhi) const
Fills the arguments with the boundaries of the world.
Float_t x1[n_points_granero]
Definition: compare.C:5
Double_t z
Definition: plot.C:276
Float_t y2[n_points_geant4]
Definition: compare.C:26
int fShowerInputs
Number of shower inputs to process from.
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
void ProjectToBoxEdge(const double xyz[], const double dxyz[], const double xlo, const double xhi, const double ylo, const double yhi, const double zlo, const double zhi, double xyzout[])
Propagates a point back to the surface of a box.
CLHEP::HepRandomEngine & fPoisEngine
double wrapvarBoxNo(const double var, const double low, const double high, int &boxno)
std::vector< int > fMaxShowers
double fToffset_corsika
Timing offset to account for propagation time through atmosphere, populated from db.
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
#define MF_LOG_DEBUG(id)
double fRandomXZShift
Each shower will be shifted by a random amount in xz so that showers won&#39;t repeatedly sample the same...
double fShowerBounds[6]
Boundaries of area over which showers are to be distributed (x(min), x(max), unused, y, z(min), z(max) )
Float_t x2[n_points_geant4]
Definition: compare.C:26
std::vector< double > fNShowersPerEvent
Number of showers to put in each event of duration fSampleTime; one per showerinput.
CLHEP::HepRandomEngine & fGenEngine
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::unique_ptr< Worker > art::ModuleBase::makeWorker ( WorkerParams const &  wp)
inherited

Definition at line 37 of file ModuleBase.cc.

References art::ModuleBase::doMakeWorker(), and art::NumBranchTypes.

38  {
39  return doMakeWorker(wp);
40  }
virtual std::unique_ptr< Worker > doMakeWorker(WorkerParams const &wp)=0
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::mayConsume ( InputTag const &  tag)
protectedinherited

Definition at line 82 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsume().

83  {
84  return collector_.mayConsume<T, BT>(tag);
85  }
ProductToken< T > mayConsume(InputTag const &)
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename T , BranchType BT>
void art::ModuleBase::mayConsumeMany ( )
protectedinherited

Definition at line 96 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeMany().

97  {
98  collector_.mayConsumeMany<T, BT>();
99  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const &  tag)
inherited

Definition at line 89 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeView().

90  {
91  return collector_.mayConsumeView<T, BT>(tag);
92  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
ViewToken< Element > mayConsumeView(InputTag const &)
ModuleDescription const & art::ModuleBase::moduleDescription ( ) const
inherited

Definition at line 13 of file ModuleBase.cc.

References art::errors::LogicError.

Referenced by art::OutputModule::doRespondToOpenInputFile(), art::OutputModule::doWriteEvent(), art::Modifier::fillProductDescriptions(), art::OutputModule::makePlugins_(), art::OutputWorker::OutputWorker(), reco::shower::LArPandoraModularShowerCreation::produce(), art::Modifier::registerProducts(), and art::OutputModule::registerProducts().

14  {
15  if (md_.has_value()) {
16  return *md_;
17  }
18 
20  "There was an error while calling moduleDescription().\n"}
21  << "The moduleDescription() base-class member function cannot be called\n"
22  "during module construction. To determine which module is "
23  "responsible\n"
24  "for calling it, find the '<module type>:<module "
25  "label>@Construction'\n"
26  "tag in the message prefix above. Please contact artists@fnal.gov\n"
27  "for guidance.\n";
28  }
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::optional< ModuleDescription > md_
Definition: ModuleBase.h:55
void evgen::CORSIKAGen::openDBs ( std::string const &  module_label)
private

Definition at line 402 of file CORSIKAGen_module.cc.

References fdb, fGenEngine, fIFDH, fShowerCopyType, fShowerInputFiles, fShowerInputs, and MF_LOG_DEBUG.

Referenced by CORSIKAGen().

403  {
404  //choose files based on fShowerInputFiles, copy them with ifdh, open them
405  // for c2: statement is unused
406  //sqlite3_stmt *statement;
407  CLHEP::RandFlat flat(fGenEngine);
408 
409  //setup ifdh object
410  if (!fIFDH) fIFDH = new ifdh_ns::ifdh;
411  const char* ifdh_debug_env = std::getenv("IFDH_DEBUG_LEVEL");
412  if (ifdh_debug_env) {
413  mf::LogInfo("CORSIKAGen") << "IFDH_DEBUG_LEVEL: " << ifdh_debug_env << "\n";
414  fIFDH->set_debug(ifdh_debug_env);
415  }
416 
417  //get ifdh path for each file in fShowerInputFiles, put into selectedflist
418  //if 1 file returned, use that file
419  //if >1 file returned, randomly select one file
420  //if 0 returned, make exeption for missing files
421  std::vector<std::pair<std::string, long>> selectedflist;
422  for (int i = 0; i < fShowerInputs; i++) {
423  if (fShowerInputFiles[i].find("*") == std::string::npos) {
424  //if there are no wildcards, don't call findMatchingFiles
425  selectedflist.push_back(std::make_pair(fShowerInputFiles[i], 0));
426  mf::LogInfo("CorsikaGen") << "Selected" << selectedflist.back().first << "\n";
427  }
428  else {
429  //use findMatchingFiles
430  //default to IFDH approach when wildcards in path
431  std::vector<std::pair<std::string, long>> flist;
432  std::string path(gSystem->DirName(fShowerInputFiles[i].c_str()));
433  std::string pattern(gSystem->BaseName(fShowerInputFiles[i].c_str()));
434  flist = fIFDH->findMatchingFiles(path, pattern);
435  unsigned int selIndex = -1;
436  if (flist.size() == 1) { //0th element is the search path:pattern
437  selIndex = 0;
438  }
439  else if (flist.size() > 1) {
440  selIndex = (unsigned int)(flat() * (flist.size() - 1) +
441  0.5); //rnd with rounding, dont allow picking the 0th element
442  }
443  else {
444  throw cet::exception("CORSIKAGen")
445  << "No files returned for path:pattern: " << path << ":" << pattern << std::endl;
446  }
447  selectedflist.push_back(flist[selIndex]);
448  mf::LogInfo("CorsikaGen") << "For " << fShowerInputFiles[i] << ":" << pattern << "\nFound "
449  << flist.size() << " candidate files"
450  << "\nChoosing file number " << selIndex << "\n"
451  << "\nSelected " << selectedflist.back().first << "\n";
452  }
453  }
454 
455  //do the fetching, store local filepaths in locallist
456  std::vector<std::string> locallist;
457  for (unsigned int i = 0; i < selectedflist.size(); i++) {
458  mf::LogInfo("CorsikaGen") << "Fetching: " << selectedflist[i].first << " "
459  << selectedflist[i].second << "\n";
460  if (fShowerCopyType == "IFDH") {
461  std::string fetchedfile(fIFDH->fetchInput(selectedflist[i].first));
462  MF_LOG_DEBUG("CorsikaGen") << "Fetched; local path: " << fetchedfile << "; copy type: IFDH";
463  locallist.push_back(fetchedfile);
464  }
465  else if (fShowerCopyType == "DIRECT") {
466  std::string fetchedfile(selectedflist[i].first);
467  MF_LOG_DEBUG("CorsikaGen")
468  << "Fetched; local path: " << fetchedfile << "; copy type: DIRECT";
469  locallist.push_back(fetchedfile);
470  }
471  else
472  throw cet::exception("CORSIKAGen")
473  << "Error copying shower db: ShowerCopyType must be \"IFDH\" or \"DIRECT\"\n";
474  }
475 
476  //open the files in fShowerInputFilesLocalPaths with sqlite3
477  for (unsigned int i = 0; i < locallist.size(); i++) {
478  //prepare and execute statement to attach db file
479  int res = sqlite3_open(locallist[i].c_str(), &fdb[i]);
480  if (res != SQLITE_OK)
481  throw cet::exception("CORSIKAGen")
482  << "Error opening db: (" << locallist[i].c_str() << ") (" << res
483  << "): " << sqlite3_errmsg(fdb[i]) << "; memory used:<<" << sqlite3_memory_used() << "/"
484  << sqlite3_memory_highwater(0) << "\n";
485  else
486  mf::LogInfo("CORSIKAGen") << "Attached db " << locallist[i] << "\n";
487  }
488  }
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< std::string > fShowerInputFiles
Set of CORSIKA shower data files to use.
ifdh_ns::ifdh * fIFDH
(optional) flux file handling
int fShowerInputs
Number of shower inputs to process from.
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
std::string fShowerCopyType
#define MF_LOG_DEBUG(id)
CLHEP::HepRandomEngine & fGenEngine
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void evgen::CORSIKAGen::populateNShowers ( )
private

Definition at line 537 of file CORSIKAGen_module.cc.

References fBuffBox, fdb, fMaxShowers, fNShowersPerEvent, fRandomXZShift, fSampleTime, fShowerAreaExtension, fShowerBounds, fShowerFluxConstants, fShowerInputs, fToffset, and geo::GeometryCore::Iterate().

Referenced by CORSIKAGen().

538  {
539  //populate vector of the number of showers per event based on:
540  //AREA the showers are being distributed over
541  //TIME of the event (fSampleTime)
542  //flux constants that determine the overall normalizations (fShowerFluxConstants)
543  //Energy range over which the sample was generated (ERANGE_*)
544  //power spectrum over which the sample was generated (ESLOPE)
545 
546  //compute shower area based on the maximal x,z dimensions of cryostat boundaries + fShowerAreaExtension
548  for (auto const& cryostat : geom->Iterate<geo::CryostatGeo>()) {
549  double bounds[6] = {0.};
550  cryostat.Boundaries(bounds);
551  for (unsigned int bnd = 0; bnd < 6; bnd++) {
552  mf::LogVerbatim("CORSIKAGen") << "Cryo Boundary: " << bnd << "=" << bounds[bnd]
553  << " ( + Buffer=" << fBuffBox[bnd] << ")\n";
554  if (fabs(bounds[bnd]) > fabs(fShowerBounds[bnd])) { fShowerBounds[bnd] = bounds[bnd]; }
555  }
556  }
557  //add on fShowerAreaExtension without being clever
562 
563  double showersArea = (fShowerBounds[1] / 100 - fShowerBounds[0] / 100) *
564  (fShowerBounds[5] / 100 - fShowerBounds[4] / 100);
565 
566  mf::LogInfo("CORSIKAGen") << "Area extended by : " << fShowerAreaExtension
567  << "\nShowers to be distributed betweeen: x=" << fShowerBounds[0]
568  << "," << fShowerBounds[1] << " & z=" << fShowerBounds[4] << ","
569  << fShowerBounds[5]
570  << "\nShowers to be distributed with random XZ shift: "
571  << fRandomXZShift << " cm"
572  << "\nShowers to be distributed over area: " << showersArea << " m^2"
573  << "\nShowers to be distributed over time: " << fSampleTime << " s"
574  << "\nShowers to be distributed with time offset: " << fToffset
575  << " s"
576  << "\nShowers to be distributed at y: " << fShowerBounds[3] << " cm";
577 
578  //db variables
579  sqlite3_stmt* statement;
580  const std::string kStatement("select erange_high,erange_low,eslope,nshow from input");
581  double upperLimitOfEnergyRange = 0., lowerLimitOfEnergyRange = 0., energySlope = 0.,
582  oneMinusGamma = 0., EiToOneMinusGamma = 0., EfToOneMinusGamma = 0.;
583 
584  for (int i = 0; i < fShowerInputs; i++) {
585  //build and do query to get run info from databases
586  // double thisrnd=flat();//need a new random number for each query
587  if (sqlite3_prepare(fdb[i], kStatement.c_str(), -1, &statement, 0) == SQLITE_OK) {
588  int res = 0;
589  res = sqlite3_step(statement);
590  if (res == SQLITE_ROW) {
591  upperLimitOfEnergyRange = sqlite3_column_double(statement, 0);
592  lowerLimitOfEnergyRange = sqlite3_column_double(statement, 1);
593  energySlope = sqlite3_column_double(statement, 2);
594  fMaxShowers.push_back(sqlite3_column_int(statement, 3));
595  oneMinusGamma = 1 + energySlope;
596  EiToOneMinusGamma = pow(lowerLimitOfEnergyRange, oneMinusGamma);
597  EfToOneMinusGamma = pow(upperLimitOfEnergyRange, oneMinusGamma);
598  mf::LogVerbatim("CORSIKAGen")
599  << "For showers input " << i << " found e_hi=" << upperLimitOfEnergyRange
600  << ", e_lo=" << lowerLimitOfEnergyRange << ", slope=" << energySlope
601  << ", k=" << fShowerFluxConstants[i] << "\n";
602  }
603  else {
604  throw cet::exception("CORSIKAGen")
605  << "Unexpected sqlite3_step return value: (" << res << "); "
606  << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
607  }
608  }
609  else {
610  throw cet::exception("CORSIKAGen") << "Error preparing statement: (" << kStatement << "); "
611  << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
612  }
613 
614  //this is computed, how?
615  double NShowers = (M_PI * showersArea * fShowerFluxConstants[i] *
616  (EfToOneMinusGamma - EiToOneMinusGamma) / oneMinusGamma) *
617  fSampleTime;
618  fNShowersPerEvent.push_back(NShowers);
619  mf::LogVerbatim("CORSIKAGen")
620  << "For showers input " << i << " the number of showers per event is " << (int)NShowers
621  << "\n";
622  }
623  }
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
double fToffset
Time offset of sample, defaults to zero (no offset) [s].
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
double fSampleTime
Duration of sample [s].
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
int fShowerInputs
Number of shower inputs to process from.
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
std::vector< int > fMaxShowers
std::vector< double > fShowerFluxConstants
Set of flux constants to be associated with each shower data file.
double fRandomXZShift
Each shower will be shifted by a random amount in xz so that showers won&#39;t repeatedly sample the same...
double fShowerBounds[6]
Boundaries of area over which showers are to be distributed (x(min), x(max), unused, y, z(min), z(max) )
std::vector< double > fNShowersPerEvent
Number of showers to put in each event of duration fSampleTime; one per showerinput.
double fShowerAreaExtension
Extend distribution of corsika particles in x,z by this much (e.g. 1000 will extend 10 m in -x...
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< double > fBuffBox
Buffer box extensions to cryostat in each direction (6 of them: x_lo,x_hi,y_lo,y_hi,z_lo,z_hi) [cm].
void evgen::CORSIKAGen::populateTOffset ( )
private

Definition at line 503 of file CORSIKAGen_module.cc.

References fdb, fShowerInputs, and fToffset_corsika.

Referenced by CORSIKAGen().

504  {
505  //populate TOffset_corsika by finding minimum ParticleTime from db file
506 
507  sqlite3_stmt* statement;
508  const std::string kStatement("select min(t) from particles");
509  double t = 0.;
510 
511  for (int i = 0; i < fShowerInputs; i++) {
512  //build and do query to get run min(t) from each db
513  if (sqlite3_prepare(fdb[i], kStatement.c_str(), -1, &statement, 0) == SQLITE_OK) {
514  int res = 0;
515  res = sqlite3_step(statement);
516  if (res == SQLITE_ROW) {
517  t = sqlite3_column_double(statement, 0);
518  mf::LogInfo("CORSIKAGen")
519  << "For showers input " << i << " found particles.min(t)=" << t << "\n";
520  if (i == 0 || t < fToffset_corsika) fToffset_corsika = t;
521  }
522  else {
523  throw cet::exception("CORSIKAGen")
524  << "Unexpected sqlite3_step return value: (" << res << "); "
525  << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
526  }
527  }
528  else {
529  throw cet::exception("CORSIKAGen") << "Error preparing statement: (" << kStatement << "); "
530  << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
531  }
532  }
533 
534  mf::LogInfo("CORSIKAGen") << "Found corsika timeoffset [ns]: " << fToffset_corsika << "\n";
535  }
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
int fShowerInputs
Number of shower inputs to process from.
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
double fToffset_corsika
Timing offset to account for propagation time through atmosphere, populated from db.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void evgen::CORSIKAGen::produce ( art::Event evt)
virtual

Implements art::EDProducer.

Definition at line 807 of file CORSIKAGen_module.cc.

References simb::MCTruth::Add(), DEFINE_ART_MODULE, fBuffBox, simb::MCTruth::GetParticle(), GetSample(), geo::GeometryCore::Iterate(), simb::kCosmicRay, simb::MCParticle::Momentum(), simb::MCTruth::NParticles(), simb::MCParticle::Position(), art::Event::put(), and simb::MCTruth::SetOrigin().

808  {
809  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
810 
812 
813  simb::MCTruth truth;
815 
816  simb::MCTruth pretruth;
817  GetSample(pretruth);
818  mf::LogInfo("CORSIKAGen") << "GetSample number of particles returned: " << pretruth.NParticles()
819  << "\n";
820  // loop over particles in the truth object
821  for (int i = 0; i < pretruth.NParticles(); ++i) {
822  simb::MCParticle particle = pretruth.GetParticle(i);
823  const TLorentzVector& v4 = particle.Position();
824  const TLorentzVector& p4 = particle.Momentum();
825  double x0[3] = {v4.X(), v4.Y(), v4.Z()};
826  double dx[3] = {p4.Px(), p4.Py(), p4.Pz()};
827 
828  // now check if the particle goes through any cryostat in the detector
829  // if so, add it to the truth object.
830  for (auto const& cryostat : geom->Iterate<geo::CryostatGeo>()) {
831  double bounds[6] = {0.};
832  cryostat.Boundaries(bounds);
833 
834  //add a buffer box around the cryostat bounds to increase the acceptance and account for scattering
835  //By default, the buffer box has zero size
836  for (unsigned int cb = 0; cb < 6; cb++)
837  bounds[cb] = bounds[cb] + fBuffBox[cb];
838 
839  //calculate the intersection point with each cryostat surface
840  bool intersects_cryo = false;
841  for (int bnd = 0; bnd != 6; ++bnd) {
842  if (bnd < 2) {
843  double p2[3] = {bounds[bnd],
844  x0[1] + (dx[1] / dx[0]) * (bounds[bnd] - x0[0]),
845  x0[2] + (dx[2] / dx[0]) * (bounds[bnd] - x0[0])};
846  if (p2[1] >= bounds[2] && p2[1] <= bounds[3] && p2[2] >= bounds[4] &&
847  p2[2] <= bounds[5]) {
848  intersects_cryo = true;
849  break;
850  }
851  }
852  else if (bnd >= 2 && bnd < 4) {
853  double p2[3] = {x0[0] + (dx[0] / dx[1]) * (bounds[bnd] - x0[1]),
854  bounds[bnd],
855  x0[2] + (dx[2] / dx[1]) * (bounds[bnd] - x0[1])};
856  if (p2[0] >= bounds[0] && p2[0] <= bounds[1] && p2[2] >= bounds[4] &&
857  p2[2] <= bounds[5]) {
858  intersects_cryo = true;
859  break;
860  }
861  }
862  else if (bnd >= 4) {
863  double p2[3] = {x0[0] + (dx[0] / dx[2]) * (bounds[bnd] - x0[2]),
864  x0[1] + (dx[1] / dx[2]) * (bounds[bnd] - x0[2]),
865  bounds[bnd]};
866  if (p2[0] >= bounds[0] && p2[0] <= bounds[1] && p2[1] >= bounds[2] &&
867  p2[1] <= bounds[3]) {
868  intersects_cryo = true;
869  break;
870  }
871  }
872  }
873 
874  if (intersects_cryo) {
875  truth.Add(particle);
876  break; //leave loop over cryostats to avoid adding particle multiple times
877  } // end if particle goes into a cryostat
878  } // end loop over cryostats in the detector
879 
880  } // loop on particles
881 
882  mf::LogInfo("CORSIKAGen")
883  << "Number of particles from getsample crossing cryostat + bounding box: "
884  << truth.NParticles() << "\n";
885 
886  truthcol->push_back(truth);
887  evt.put(std::move(truthcol));
888 
889  return;
890  } // end produce
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:220
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
int NParticles() const
Definition: MCTruth.h:75
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:221
Event generator information.
Definition: MCTruth.h:32
void GetSample(simb::MCTruth &)
Cosmic rays.
Definition: MCTruth.h:24
std::vector< double > fBuffBox
Buffer box extensions to cryostat in each direction (6 of them: x_lo,x_hi,y_lo,y_hi,z_lo,z_hi) [cm].
void evgen::CORSIKAGen::ProjectToBoxEdge ( const double  xyz[],
const double  dxyz[],
const double  xlo,
const double  xhi,
const double  ylo,
const double  yhi,
const double  zlo,
const double  zhi,
double  xyzout[] 
)
private

Propagates a point back to the surface of a box.

Parameters
xyzcoordinates of the point to be propagated ({ x, y, z })
dxyzdirection of the point ({ dx, dy, dz })
xlolower x coordinate of the target box
xhiupper x coordinate of the target box
ylolower y coordinate of the target box
yhiupper y coordinate of the target box
zlolower z coordinate of the target box
zhiupper z coordinate of the target box
xyzout_(output, room for at least 3 numbers)_ propagated point

The point xyz, assumed to be inside the box, is propagated at the level of the closest among the sides of the box. Note that this means the propagated point might still be not on the surface of the box, even if it would later reach it. It is anyway guaranteed that xyzout is not inside the target box, and that at least one of its three coordinates is on the box surface.

Definition at line 356 of file CORSIKAGen_module.cc.

References d.

Referenced by GetSample().

365  {
366 
367  //we want to project backwards, so take mirror of momentum
368  const double dxyz[3] = {-indxyz[0], -indxyz[1], -indxyz[2]};
369 
370  // Compute the distances to the x/y/z walls
371  double dx = 99.E99;
372  double dy = 99.E99;
373  double dz = 99.E99;
374  if (dxyz[0] > 0.0) { dx = (xhi - xyz[0]) / dxyz[0]; }
375  else if (dxyz[0] < 0.0) {
376  dx = (xlo - xyz[0]) / dxyz[0];
377  }
378  if (dxyz[1] > 0.0) { dy = (yhi - xyz[1]) / dxyz[1]; }
379  else if (dxyz[1] < 0.0) {
380  dy = (ylo - xyz[1]) / dxyz[1];
381  }
382  if (dxyz[2] > 0.0) { dz = (zhi - xyz[2]) / dxyz[2]; }
383  else if (dxyz[2] < 0.0) {
384  dz = (zlo - xyz[2]) / dxyz[2];
385  }
386 
387  // Choose the shortest distance
388  double d = 0.0;
389  if (dx < dy && dx < dz)
390  d = dx;
391  else if (dy < dz && dy < dx)
392  d = dy;
393  else if (dz < dx && dz < dy)
394  d = dz;
395 
396  // Make the step
397  for (int i = 0; i < 3; ++i) {
398  xyzout[i] = xyz[i] + dxyz[i] * d;
399  }
400  }
Float_t d
Definition: plot.C:235
void art::Modifier::registerProducts ( ProductDescriptions productsToRegister)
inherited

Definition at line 16 of file Modifier.cc.

References art::ModuleBase::moduleDescription(), and art::ProductRegistryHelper::registerProducts().

17  {
18  ProductRegistryHelper::registerProducts(productsToRegister,
20  }
void registerProducts(ProductDescriptions &productsToRegister, ModuleDescription const &md)
ModuleDescription const & moduleDescription() const
Definition: ModuleBase.cc:13
void art::ModuleBase::setModuleDescription ( ModuleDescription const &  md)
inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

32  {
33  md_ = md;
34  }
std::optional< ModuleDescription > md_
Definition: ModuleBase.h:55
void art::ModuleBase::sortConsumables ( std::string const &  current_process_name)
inherited

Definition at line 49 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::sortConsumables().

50  {
51  // Now that we know we have seen all the consumes declarations,
52  // sort the results for fast lookup later.
53  collector_.sortConsumables(current_process_name);
54  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
void sortConsumables(std::string const &current_process_name)
double evgen::CORSIKAGen::wrapvar ( const double  var,
const double  low,
const double  high 
)
private

Definition at line 490 of file CORSIKAGen_module.cc.

Referenced by GetSample().

491  {
492  //wrap variable so that it's always between low and high
493  return (var - (high - low) * floor(var / (high - low))) + low;
494  }
double evgen::CORSIKAGen::wrapvarBoxNo ( const double  var,
const double  low,
const double  high,
int &  boxno 
)
private

Definition at line 496 of file CORSIKAGen_module.cc.

Referenced by GetSample().

497  {
498  //wrap variable so that it's always between low and high
499  boxno = int(floor(var / (high - low)));
500  return (var - (high - low) * floor(var / (high - low))) + low;
501  }

Member Data Documentation

std::vector<double> evgen::CORSIKAGen::fBuffBox
private

Buffer box extensions to cryostat in each direction (6 of them: x_lo,x_hi,y_lo,y_hi,z_lo,z_hi) [cm].

Definition at line 289 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), populateNShowers(), and produce().

sqlite3* evgen::CORSIKAGen::fdb[5]
private

Pointers to sqlite3 database object, max of 5.

Definition at line 292 of file CORSIKAGen_module.cc.

Referenced by GetSample(), openDBs(), populateNShowers(), populateTOffset(), and ~CORSIKAGen().

CLHEP::HepRandomEngine& evgen::CORSIKAGen::fGenEngine
private

Definition at line 295 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), GetSample(), and openDBs().

ifdh_ns::ifdh* evgen::CORSIKAGen::fIFDH = 0
private

(optional) flux file handling

Definition at line 278 of file CORSIKAGen_module.cc.

Referenced by openDBs(), and ~CORSIKAGen().

std::vector<int> evgen::CORSIKAGen::fMaxShowers
private

Definition at line 268 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

std::vector<double> evgen::CORSIKAGen::fNShowersPerEvent
private

Number of showers to put in each event of duration fSampleTime; one per showerinput.

Definition at line 267 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

CLHEP::HepRandomEngine& evgen::CORSIKAGen::fPoisEngine
private

Definition at line 296 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), and GetSample().

double evgen::CORSIKAGen::fProjectToHeight = 0.
private

Height to which particles will be projected [cm].

Definition at line 281 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), and GetSample().

double evgen::CORSIKAGen::fRandomXZShift
private
Initial value:
=
0.

Each shower will be shifted by a random amount in xz so that showers won't repeatedly sample the same space [cm].

Definition at line 293 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), GetSample(), and populateNShowers().

double evgen::CORSIKAGen::fSampleTime = 0.
private

Duration of sample [s].

Definition at line 286 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), GetSample(), and populateNShowers().

double evgen::CORSIKAGen::fShowerAreaExtension
private
Initial value:
=
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) [cm].

Definition at line 290 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), and populateNShowers().

double evgen::CORSIKAGen::fShowerBounds[6]
private
Initial value:
= {
0.,
0.,
0.,
0.,
0.,
0.}

Boundaries of area over which showers are to be distributed (x(min), x(max), unused, y, z(min), z(max) )

Definition at line 269 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

std::string evgen::CORSIKAGen::fShowerCopyType
private

Definition at line 283 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), and openDBs().

std::vector<double> evgen::CORSIKAGen::fShowerFluxConstants
private

Set of flux constants to be associated with each shower data file.

Definition at line 285 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), and populateNShowers().

std::vector<std::string> evgen::CORSIKAGen::fShowerInputFiles
private

Set of CORSIKA shower data files to use.

Definition at line 282 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), and openDBs().

int evgen::CORSIKAGen::fShowerInputs = 0
private

Number of shower inputs to process from.

Definition at line 265 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), GetSample(), openDBs(), populateNShowers(), populateTOffset(), and ~CORSIKAGen().

double evgen::CORSIKAGen::fToffset = 0.
private

Time offset of sample, defaults to zero (no offset) [s].

Definition at line 287 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), GetSample(), and populateNShowers().

double evgen::CORSIKAGen::fToffset_corsika
private
Initial value:
=
0.

Timing offset to account for propagation time through atmosphere, populated from db.

Definition at line 276 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateTOffset().


The documentation for this class was generated from the following file: