Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes | List of all members
evgen::CORSIKAGen Class Reference
Inheritance diagram for evgen::CORSIKAGen:
art::EDProducer art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 CORSIKAGen (fhicl::ParameterSet const &p)
 
virtual ~CORSIKAGen ()
 
void produce (art::Event &evt) override
 
virtual void beginRun (art::Run &run)
 
virtual void beginSubRun (art::SubRun &run)
 
virtual void endSubRun (art::SubRun &run)
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Private Member Functions

void openDBs ()
 
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 DetectorBigBoxCut (simb::MCTruth &truth)
 
bool isIntersectTheBox (const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
 
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[])
 

Private Attributes

TStopwatch stopwatch_
 keep track of how long it takes to run the job More...
 
genie::RandomGenrnd_
 Random generator from GENIE. More...
 
int m_ShowerInputs = 0
 Number of shower inputs to process from. More...
 
std::vector< double > m_NShowersPerEvent
 duration m_fcl_SampleTime; one per showerinput More...
 
std::vector< intm_MaxShowers
 
double m_ShowerBounds [6]
 
double m_Toffset_corsika = 0.
 atmosphere, populated from db (optional) flux file handling More...
 
ifdh_ns::ifdh * m_IFDH = 0
 (optional) flux file handling More...
 
sqlite3 * fdb [5]
 Pointers to sqlite3 database object, max of 5. More...
 
double m_fcl_ProjectToHeight = 0.
 Height to which particles will be projected [cm]. More...
 
std::vector< std::stringm_fcl_ShowerInputFiles
 Set of CORSIKA shower data files to use. More...
 
std::vector< double > m_fcl_ShowerFluxConstants
 with each shower data file More...
 
double m_fcl_SampleTime = 0.
 Duration of sample [s]. More...
 
double m_fcl_Toffset = 0.
 Time offset of sample, defaults to zero (no offset) [s]. More...
 
double m_fcl_ShowerAreaExtension = 0.
 (e.g. 1000 will extend 10 m in -x, +x, -z, and +z) [cm] More...
 
double m_fcl_RandomXZShift = 0.
 showers won't repeatedly sample the same space [cm] More...
 
bool m_fcl_IsBigBoxUsed
 
bool m_fcl_EnhanceNNbarBkgd
 Generating only cosmogenic (anti)neutrons and gammas if TRUE. More...
 
double m_fcl_NNbarBkgdEnergyCut
 Energy (in GeV) cut for NNbar background particles (neutron, antineutron, gamma) More...
 
int m_fcl_Cycle
 MC cycle generation number. More...
 

Detailed Description

Definition at line 54 of file CORSIKAGen_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

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

Definition at line 121 of file CORSIKAGen_module.cc.

References genie::RandomGen::Instance(), LOG_INFO, m_fcl_ProjectToHeight, m_fcl_SampleTime, m_fcl_ShowerFluxConstants, m_fcl_ShowerInputFiles, m_ShowerInputs, openDBs(), populateNShowers(), populateTOffset(), rnd_, genie::RandomGen::SetSeed(), and stopwatch_.

122  : m_fcl_ProjectToHeight(p.get<double>("ProjectToHeight", 0.)),
123  m_fcl_ShowerInputFiles(p.get<std::vector<std::string>>("ShowerInputFiles")),
124  m_fcl_ShowerFluxConstants(p.get<std::vector<double>>("ShowerFluxConstants")),
125  m_fcl_SampleTime(p.get<double>("SampleTime", 0.)),
126  m_fcl_Toffset(p.get<double>("TimeOffset", 0.)),
127  m_fcl_ShowerAreaExtension(p.get<double>("ShowerAreaExtension", 0.)),
128  m_fcl_RandomXZShift(p.get<double>("RandomXZShift", 0.)),
129  m_fcl_IsBigBoxUsed(p.get<bool>("IsBigBoxUsed", true)),
130  m_fcl_Cycle(p.get<int>("Cycle", 0)),
131  m_fcl_EnhanceNNbarBkgd(p.get<bool>("EnhanceNNbarBkgd", true)),
132  m_fcl_NNbarBkgdEnergyCut(p.get<double>("NNbarBkgdEnergyCut", 1.0)) {
133  stopwatch_.Start();
134 
136  rnd_->SetSeed(std::time(NULL));
137 
138  if (m_fcl_ShowerInputFiles.size() != m_fcl_ShowerFluxConstants.size() || m_fcl_ShowerInputFiles.size() == 0 || m_fcl_ShowerFluxConstants.size() == 0) {
139  throw cet::exception("CORSIKAGen") << "ShowerInputFiles and ShowerFluxConstants have different or invalid sizes!" << "\n";
140  }
142  if (m_fcl_SampleTime == 0.) throw cet::exception("CORSIKAGen") << "SampleTime not set!";
143  if (m_fcl_ProjectToHeight == 0.) LOG_INFO("CORSIKAGen") << "Using 0. for m_fcl_ProjectToHeight!";
144 
145  // create a default random engine; obtain the random seed from
146  // NuRandomService, unless overridden in configuration with key "Seed"
147  this->openDBs();
148  this->populateNShowers();
149  this->populateTOffset();
150  produces<std::vector<simb::MCTruth>>();
151  produces<sumdata::SubRunData, art::InSubRun>();
152  produces<sumdata::RunData, art::InRun>();
153 }
bool m_fcl_EnhanceNNbarBkgd
Generating only cosmogenic (anti)neutrons and gammas if TRUE.
double m_fcl_ProjectToHeight
Height to which particles will be projected [cm].
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
double m_fcl_RandomXZShift
showers won&#39;t repeatedly sample the same space [cm]
const char * p
Definition: xmltok.h:285
genie::RandomGen * rnd_
Random generator from GENIE.
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::vector< std::string > m_fcl_ShowerInputFiles
Set of CORSIKA shower data files to use.
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
TStopwatch stopwatch_
keep track of how long it takes to run the job
double m_fcl_Toffset
Time offset of sample, defaults to zero (no offset) [s].
int m_ShowerInputs
Number of shower inputs to process from.
double m_fcl_SampleTime
Duration of sample [s].
double m_fcl_NNbarBkgdEnergyCut
Energy (in GeV) cut for NNbar background particles (neutron, antineutron, gamma)
std::vector< double > m_fcl_ShowerFluxConstants
with each shower data file
#define LOG_INFO(stream)
Definition: Messenger.h:144
double m_fcl_ShowerAreaExtension
(e.g. 1000 will extend 10 m in -x, +x, -z, and +z) [cm]
int m_fcl_Cycle
MC cycle generation number.
void SetSeed(long int seed)
Definition: RandomGen.cxx:90
evgen::CORSIKAGen::~CORSIKAGen ( )
virtual

Definition at line 196 of file CORSIKAGen_module.cc.

References fdb, MECModelEnuComparisons::i, LOG_INFO, m_IFDH, m_ShowerInputs, and stopwatch_.

196  {
197  for (int i = 0; i < m_ShowerInputs; i++) {
198  sqlite3_close(fdb[i]);
199  }
200  // cleanup temp files
201  m_IFDH->cleanup();
202 
203  stopwatch_.Stop();
204  LOG_INFO("CORSIKAGen")
205  << "real time to produce file: " << stopwatch_.RealTime();
206 }
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
TStopwatch stopwatch_
keep track of how long it takes to run the job
int m_ShowerInputs
Number of shower inputs to process from.
ifdh_ns::ifdh * m_IFDH
(optional) flux file handling
#define LOG_INFO(stream)
Definition: Messenger.h:144

Member Function Documentation

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

Reimplemented from art::EDProducer.

Definition at line 208 of file CORSIKAGen_module.cc.

References geo::GeometryBase::DetId(), geo::GeometryBase::ExtractGDML(), geo::GeometryBase::FileBaseName(), and art::Run::put().

208  {
210  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(
211  geo->DetId(), geo->FileBaseName(), geo->ExtractGDML()));
212 
213  run.put(std::move(runcol));
214 }
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
Helper for AttenCurve.
Definition: Path.h:10
std::string FileBaseName() const
void evgen::CORSIKAGen::beginSubRun ( art::SubRun run)
virtual

Reimplemented from art::EDProducer.

Definition at line 216 of file CORSIKAGen_module.cc.

216  {
217 }
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDProducer::currentContext ( ) const
protectedinherited
void evgen::CORSIKAGen::DetectorBigBoxCut ( simb::MCTruth truth)
private

Definition at line 657 of file CORSIKAGen_module.cc.

References simb::MCTruth::Add(), geo::GeometryBase::DetectorBigBox(), geo::GeometryBase::DetId(), simb::MCParticle::EndX(), simb::MCParticle::EndY(), simb::MCParticle::EndZ(), simb::MCTruth::GetParticle(), ipart, isIntersectTheBox(), novadaq::cnv::kNEARDET, simb::MCParticle::Mass(), simb::MCParticle::Momentum(), simb::MCParticle::Mother(), npart, simb::MCTruth::NParticles(), simb::MCParticle::P(), part, simb::MCParticle::PdgCode(), simb::MCParticle::Process(), simb::MCParticle::Px(), simb::MCParticle::Py(), simb::MCParticle::Pz(), simb::MCParticle::StatusCode(), simb::MCParticle::T(), simb::MCParticle::TrackId(), x1, submit_syst::x2, y1, and submit_syst::y2.

Referenced by produce().

657  {
658  simb::MCTruth mctruth_new;
659 
660  double x1 = 0;
661  double x2 = 0;
662  double y1 = 0;
663  double y2 = 0;
664  double z1 = 0;
665  double z2 = 0;
666 
668  // Get the ranges of the x, y, and z for the "Detector Volume" that the entire
669  // detecotr geometry lives in. If any of the pointers is NULL, the
670  // corresponding coordinate is ignored.
671  geo->DetectorBigBox(&x1, &x2, &y1, &y2, &z1, &z2);
672 
673  // Check for any missing range
674  if (x1 == 0 && x2 == 0 && y1 == 0 && y2 == 0 && z1 == 0 && z2 == 0) {
675  throw cet::exception("NoGeometryBoxCoords") << "No Geometry Box Coordinates Set\n" << __FILE__ << ":" << __LINE__ << "\n";
676  return;
677  }
678 
679  // Loop over the particles stored by CRY
680  int npart = truth.NParticles();
681  for (int ipart = 0; ipart < npart; ++ipart) {
682  simb::MCParticle particle = truth.GetParticle(ipart);
683 
684  double px = particle.Px() / particle.P();
685  double py = particle.Py() / particle.P();
686  double pz = particle.Pz() / particle.P();
687 
688  double vx = particle.EndX();
689  double vy = particle.EndY();
690  double vz = particle.EndZ();
691 
692  // For the near detector, move all the cosmics from their position in
693  // the window, straight down, to 1cm above the detector big box. This
694  // ensures that a much larger fraction of them pass the IntersectTheBox
695  // cut below. Subsequently the ProjectCosmicsToSurface call will
696  // backtrack them up to the surface. The flux remains correct throughout
697  // this scheme due to the homogeneity of the flux in space and time.
698  // It does get the correlations wrong, but the chances of two correlated
699  // events making it to the ND are remote anyway.
700  if (geo->DetId() == novadaq::cnv::kNEARDET) vy = y2 + 1;
701 
702  double xyz[3] = {vx, vy, vz};
703  double dxyz[3] = {px, py, pz};
704 
705  // If intersecting the DetectorBigBox, add particle
706  // currently, using the methods in CosmicsGen code until the geometry methods can be validated
707  if (this->isIntersectTheBox(xyz, dxyz, x1, x2, y1, y2, z1, z2)) {
708  simb::MCParticle part(particle.TrackId(), particle.PdgCode(), particle.Process(), particle.Mother(), particle.Mass(), particle.StatusCode());
709  part.AddTrajectoryPoint(TLorentzVector(vx, vy, vz, particle.T()), particle.Momentum());
710  mctruth_new.Add(part);
711  }
712  } // end of loop over particles
713 
714  // reassign the MCTruth
715  truth = mctruth_new;
716 }
int PdgCode() const
Definition: MCParticle.h:211
double Py(const int i=0) const
Definition: MCParticle.h:230
double EndZ() const
Definition: MCParticle.h:227
Int_t ipart
Definition: Style.C:10
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
int Mother() const
Definition: MCParticle.h:212
double Mass() const
Definition: MCParticle.h:238
double Px(const int i=0) const
Definition: MCParticle.h:229
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
int StatusCode() const
Definition: MCParticle.h:210
int NParticles() const
Definition: MCTruth.h:74
std::string Process() const
Definition: MCParticle.h:214
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
double EndY() const
Definition: MCParticle.h:226
int TrackId() const
Definition: MCParticle.h:209
TString part[npart]
Definition: Style.C:32
Int_t npart
Definition: Style.C:7
double P(const int i=0) const
Definition: MCParticle.h:233
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
double T(const int i=0) const
Definition: MCParticle.h:223
Near Detector in the NuMI cavern.
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
double Pz(const int i=0) const
Definition: MCParticle.h:231
bool isIntersectTheBox(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
Event generator information.
Definition: MCTruth.h:32
double EndX() const
Definition: MCParticle.h:225
Helper for AttenCurve.
Definition: Path.h:10
void DetectorBigBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
void evgen::CORSIKAGen::endSubRun ( art::SubRun run)
virtual

Reimplemented from art::EDProducer.

Definition at line 219 of file CORSIKAGen_module.cc.

References m_fcl_Cycle, art::SubRun::put(), and sd().

219  {
220  // store the cycle information
221  std::unique_ptr<sumdata::SubRunData> sd(new sumdata::SubRunData(m_fcl_Cycle));
222  subrun.put(std::move(sd));
223 }
double sd(Eigen::VectorXd x)
int m_fcl_Cycle
MC cycle generation number.
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

Referenced by skim::NueSkimmer::CopyMichelSlice(), and skim::NueSkimmer::CopyMichelTrack().

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
void evgen::CORSIKAGen::GetSample ( simb::MCTruth mctruth)
private

Definition at line 435 of file CORSIKAGen_module.cc.

References simb::MCTruth::Add(), simb::MCParticle::AddTrajectoryPoint(), dx, fdb, genie::utils::style::Format(), geom(), MECModelEnuComparisons::i, LOG_INFO, m, m_fcl_ProjectToHeight, m_fcl_RandomXZShift, m_fcl_SampleTime, m_fcl_Toffset, m_MaxShowers, m_NShowersPerEvent, m_ShowerBounds, m_ShowerInputs, m_Toffset_corsika, make_root_from_grid_output::pdg, elec2geo::pos, ProjectToBoxEdge(), rnd_, genie::RandomGen::RndNum(), confusionMatrixTree::t, geo::GeometryBase::WorldBox(), wrapvar(), wrapvarBoxNo(), submit_syst::x, x1, submit_syst::x2, y1, submit_syst::y2, and test::z.

Referenced by produce().

435  {
436  // for each input, randomly pull m_NShowersPerEvent[i] showers from the
437  // Particles table and randomly place them in time (between -m_fcl_SampleTime/2
438  // and m_fcl_SampleTime/2) wrap their positions based on the size of the area
439  // under consideration based on
440  // http://nusoft.fnal.gov/larsoft/doxsvn/html/CRYHelper_8cxx_source.html
441  // (Sample) query from sqlite db with select * from particles where shower in
442  // (select id from showers ORDER BY substr(id*0.51123124141,length(id)+2) limit
443  // 100000) ORDER BY substr(shower*0.51123124141,length(shower)+2); where
444  // 0.51123124141 is a random seed to allow randomly selecting rows and should
445  // be randomly generated for each query the inner order by is to select
446  // randomly from the possible shower id's the outer order by is to make sure
447  // the shower numbers are ordered randomly (without this, the showers always
448  // come out ordered by shower number and 100000 is the number of showers to be
449  // selected at random and needs to be less than the number of showers in the
450  // showers table TDatabasePDG is for looking up particle masses
451  static TDatabasePDG *pdgt = TDatabasePDG::Instance();
452 
453  // db variables
454  sqlite3_stmt *statement;
455  const TString kStatement(
456  "select shower,pdg,px,py,pz,x,z,t,e from particles where shower in "
457  "(select id from showers ORDER BY substr(id*%f,length(id)+2) limit %d) "
458  "ORDER BY substr(shower*%f,length(shower)+2)");
459  // get geometry and figure where to project particles to, based on CRYHelper
461  double x1 = 0.;
462  double x2 = 0.;
463  double y1 = 0.;
464  double y2 = 0.;
465  double z1 = 0.;
466  double z2 = 0.;
467  geom->WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);
468  // make the world box slightly smaller so that the projection to the edge avoids possible rounding errors later on with Geant4
469 
470  double fBoxDelta = 1.e-5;
471 
472  x1 += fBoxDelta;
473  x2 -= fBoxDelta;
474  y1 += fBoxDelta;
476  z1 += fBoxDelta;
477  z2 -= fBoxDelta;
478 
479  // populate mctruth
480  int ntotalCtr = 0; // count number of particles added to mctruth
481  int lastShower = 0; // keep track of last shower id so that t can be randomized on every new shower
482  long int nShowerCntr = 0; // keep track of how many showers are left to be added to mctruth
483  int nShowerQry = 0; // number of showers to query from db
484  int shower, pdg;
485  double px, py, pz, x, z, tParticleTime, etot, showerTime = 0., showerTimex = 0., showerTimez = 0., showerXOffset = 0., showerZOffset = 0., t;
486  for (int i = 0; i < m_ShowerInputs; i++) {
487  nShowerCntr = rnd_->RndNum().Poisson(m_NShowersPerEvent[i]);
488  LOG_INFO("CORSIKAGEN") << " Shower input " << i << " with mean " << m_NShowersPerEvent[i] << " generating " << nShowerCntr
489  << ". Maximum number of showers in database: " << m_MaxShowers[i];
490  while (nShowerCntr > 0) {
491  // how many showers should we query?
492  if (nShowerCntr > m_MaxShowers[i]) nShowerQry = m_MaxShowers[i]; // take the group size
493  else nShowerQry = nShowerCntr; // take the rest that are needed
494 
495  // build and do query to get nshowers
496  double thisrnd = rnd_->RndNum().Rndm(); // need a new random number for each query
497  TString kthisStatement = TString::Format(kStatement.Data(), thisrnd, nShowerQry, thisrnd);
498  LOG_INFO("CORSIKAGen") << "Executing: " << kthisStatement;
499  if (sqlite3_prepare(fdb[i], kthisStatement.Data(), -1, &statement, 0) == SQLITE_OK) {
500  int res = 0;
501  // loop over database rows, pushing particles into mctruth object
502  while (1) {
503  res = sqlite3_step(statement);
504  if (res == SQLITE_ROW) {
505  /*
506  * Memo columns:
507  * [0] shower
508  * [1] particle ID (PDG)
509  * [2] momentum: x component [GeV/c]
510  * [3] momentum: y component [GeV/c]
511  * [4] momentum: z component [GeV/c]
512  * [5] position: x component [cm]
513  * [6] position: z component [cm]
514  * [7] time [ns]
515  * [8] energy [GeV]
516  */
517  shower = sqlite3_column_int(statement, 0);
518  if (shower != lastShower) {
519  // each new shower gets its own random time and position offsets
520  showerTime = 1e9 * (rnd_->RndNum().Rndm() * m_fcl_SampleTime); // converting from s to ns
521  showerTimex = 1e9 * (rnd_->RndNum().Rndm() * m_fcl_SampleTime); // converting from s to ns
522  showerTimez = 1e9 * (rnd_->RndNum().Rndm() * m_fcl_SampleTime); // converting from s to ns
523  // and a random offset in both z and x controlled by the
524  // m_fcl_RandomXZShift parameter
525  showerXOffset = rnd_->RndNum().Rndm() * m_fcl_RandomXZShift - (m_fcl_RandomXZShift / 2);
526  showerZOffset = rnd_->RndNum().Rndm() * m_fcl_RandomXZShift - (m_fcl_RandomXZShift / 2);
527  }
528  pdg = sqlite3_column_int(statement, 1);
529  // get mass for this particle
530  double m = 0.; // in GeV
531  TParticlePDG *pdgp = pdgt->GetParticle(pdg);
532  if (pdgp) m = pdgp->Mass();
533  // Note: position/momentum in db have north=-x and west=+z, rotate
534  // so that +z is north and +x is west get momentum components
535  // To do: Can we use the same orientation for NOvA? My uneducated guess is yes. Check with detsim.
536  px = sqlite3_column_double(statement, 4); // uboone x=Particlez
537  py = sqlite3_column_double(statement, 3);
538  pz = -sqlite3_column_double(statement, 2); // uboone z=-Particlex
539  etot = sqlite3_column_double(statement, 8);
540 
541  // get/calculate position components
542  int boxnoX = 0, boxnoZ = 0;
543  x = wrapvarBoxNo(sqlite3_column_double(statement, 6) + showerXOffset, m_ShowerBounds[0], m_ShowerBounds[1], boxnoX);
544  z = wrapvarBoxNo(-sqlite3_column_double(statement, 5) + showerZOffset, m_ShowerBounds[4], m_ShowerBounds[5], boxnoZ);
545  tParticleTime = sqlite3_column_double(statement, 7); // time offset, includes propagation time from top of atmosphere
546  // actual particle time is particle surface arrival time
547  //+ shower start time
548  //+ global offset (fcl parameter, in s)
549  //- propagation time through atmosphere
550  //+ boxNo{X,Z} time offset to make grid boxes have different shower times
551  t = tParticleTime + showerTime + (1e9 * m_fcl_Toffset) - m_Toffset_corsika + showerTimex * boxnoX + showerTimez * boxnoZ;
552  // wrap surface arrival so that it's in the desired time window
553  t = wrapvar(t, (1e9 * m_fcl_Toffset), 1e9 * (m_fcl_Toffset + m_fcl_SampleTime));
554  simb::MCParticle p(ntotalCtr, pdg, "primary", -200, m, 1);
555 
556  // project back to worldvol/m_fcl_ProjectToHeight
557  /*
558  * This back propagation goes from a point on the upper surface of
559  * the cryostat back to the edge of the world, except that that
560  * world is cut short by `m_fcl_ProjectToHeight` (`y2`) ceiling.
561  * The projection will most often lie on that ceiling, but it may
562  * end up instead on one of the side edges of the world, or even
563  * outside it.
564  */
565  double xyzo[3];
566  double x0[3] = {x, m_ShowerBounds[3], z};
567  double dx[3] = {px, py, pz};
568  this->ProjectToBoxEdge(x0, dx, x1, x2, y1, y2, z1, z2, xyzo);
569  TLorentzVector pos(xyzo[0], xyzo[1], xyzo[2], t); // time needs to be in ns to match GENIE, etc
570  TLorentzVector mom(px, py, pz, etot);
571  p.AddTrajectoryPoint(pos, mom);
572  mctruth.Add(p);
573  ntotalCtr++;
574  lastShower = shower;
575  } else if (res == SQLITE_DONE) {
576  break;
577  } else {
578  throw cet::exception("CORSIKAGen") << "Unexpected sqlite3_step return value: (" << res << "); " << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
579  }
580  } // end while loop particle query
581  } else {
582  throw cet::exception("CORSIKAGen") << "Error preparing statement: (" << kthisStatement << "); " << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
583  }
584  nShowerCntr = nShowerCntr - nShowerQry;
585  } // end while loop counting nShowerCntr
586  }
587 }
double wrapvar(const double var, const double low, const double high)
double m_fcl_ProjectToHeight
Height to which particles will be projected [cm].
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
double m_fcl_RandomXZShift
showers won&#39;t repeatedly sample the same space [cm]
std::vector< int > m_MaxShowers
const char * p
Definition: xmltok.h:285
genie::RandomGen * rnd_
Random generator from GENIE.
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
double dx[NP][NC]
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
void WorldBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
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[])
std::vector< double > m_NShowersPerEvent
duration m_fcl_SampleTime; one per showerinput
double m_Toffset_corsika
atmosphere, populated from db (optional) flux file handling
double m_fcl_Toffset
Time offset of sample, defaults to zero (no offset) [s].
int m_ShowerInputs
Number of shower inputs to process from.
double wrapvarBoxNo(const double var, const double low, const double high, int &boxno)
z
Definition: test.py:28
TRandom3 & RndNum(void) const
rnd number generator used by MC integrators & other numerical methods
Definition: RandomGen.h:78
double m_fcl_SampleTime
Duration of sample [s].
void geom(int which=0)
Definition: geom.C:163
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
#define LOG_INFO(stream)
Definition: Messenger.h:144
bool evgen::CORSIKAGen::isIntersectTheBox ( const double  xyz[],
const double  dxyz[],
double  xlo,
double  xhi,
double  ylo,
double  yhi,
double  zlo,
double  zhi 
)
private

Definition at line 718 of file CORSIKAGen_module.cc.

References DEFINE_ART_MODULE(), dx, dy, dz, stan::math::fabs(), submit_syst::x, xhi, make_syst_table_plots::xlo, submit_syst::y, yhi, ylo, and test::z.

Referenced by DetectorBigBoxCut().

721  {
722  // If inside the box, obviously intesected
723  if (xyz[0] >= xlo && xyz[0] <= xhi && xyz[1] >= ylo && xyz[1] <= yhi && xyz[2] >= zlo && xyz[2] <= zhi)
724  return true;
725 
726  // So, now we know that the particle is outside the box
727  double dx = 0., dy = 0., dz = 0.;
728 
729  // Checking intersection with 6 planes
730  // 1. Check intersection with the upper plane
731  dy = xyz[1] - yhi;
732  // Check whether the track going from above and down or from below and up.
733  // Otherwise not going to intersect this plane
734  if ((dy > 0 && dxyz[1] < 0) || (dy < 0 && dxyz[1] > 0)) {
735  double dl = fabs(dy / dxyz[1]);
736  double x = xyz[0] + dxyz[0] * dl;
737  double z = xyz[2] + dxyz[2] * dl;
738 
739  // is it inside the square?
740  if (x >= xlo && x <= xhi && z >= zlo && z <= zhi) return true;
741  }
742 
743  // 2. Check intersection with the lower plane
744  dy = xyz[1] - ylo;
745  // Check whether the track going from above and down or from below and up.
746  // Otherwise not going to intersect this plane
747  if ((dy > 0 && dxyz[1] < 0) || (dy < 0 && dxyz[1] > 0)) {
748  double dl = fabs(dy / dxyz[1]);
749  double x = xyz[0] + dxyz[0] * dl;
750  double z = xyz[2] + dxyz[2] * dl;
751 
752  // is it inside the square?
753  if (x >= xlo && x <= xhi && z >= zlo && z <= zhi) return true;
754  }
755 
756  // 3. Check intersection with the right plane
757  dz = xyz[2] - zhi;
758  // Check whether the track going from right part to the left or from left part
759  // to right. Otherwise not going to intersect this plane
760  if ((dz > 0 && dxyz[2] < 0) || (dz < 0 && dxyz[2] > 0)) {
761  double dl = fabs(dz / dxyz[2]);
762  double x = xyz[0] + dxyz[0] * dl;
763  double y = xyz[1] + dxyz[1] * dl;
764 
765  // is it inside the square?
766  if (x >= xlo && x <= xhi && y >= ylo && y <= yhi) return true;
767  }
768 
769  // 4. Check intersection with the left plane
770  dz = xyz[2] - zlo;
771  // Check whether the track going from right part to the left or from left part
772  // to right. Otherwise not going to intersect this plane
773  if ((dz > 0 && dxyz[2] < 0) || (dz < 0 && dxyz[2] > 0)) {
774  double dl = fabs(dz / dxyz[2]);
775  double x = xyz[0] + dxyz[0] * dl;
776  double y = xyz[1] + dxyz[1] * dl;
777 
778  // is it inside the square?
779  if (x >= xlo && x <= xhi && y >= ylo && y <= yhi) return true;
780  }
781 
782  // 5. Check intersection with the far plane
783  dx = xyz[0] - xhi;
784  // Check whether the track going from far part toward us or from near part to
785  // right. Otherwise not going to intersect this plane
786  if ((dx > 0 && dxyz[0] < 0) || (dx < 0 && dxyz[0] > 0)) {
787  double dl = fabs(dx / dxyz[0]);
788  double y = xyz[1] + dxyz[1] * dl;
789  double z = xyz[2] + dxyz[2] * dl;
790 
791  // is it inside the square?
792  if (z >= zlo && z <= zhi && y >= ylo && y <= yhi) return true;
793  }
794 
795  // 6. Check intersection with the near plane
796  dx = xyz[0] - xlo;
797  // Check whether the track going from far part toward us or from near part to
798  // right. Otherwise not going to intersect this plane
799  if ((dx > 0 && dxyz[0] < 0) || (dx < 0 && dxyz[0] > 0)) {
800  double dl = fabs(dx / dxyz[0]);
801  double y = xyz[1] + dxyz[1] * dl;
802  double z = xyz[2] + dxyz[2] * dl;
803 
804  // is it inside the square?
805  if (z >= zlo && z <= zhi && y >= ylo && y <= yhi) return true;
806  }
807 
808  return false;
809 }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double dy[NP][NC]
double dx[NP][NC]
double dz[NP][NC]
z
Definition: test.py:28
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID(), and string.

41  {
42  return true;
43  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited
void evgen::CORSIKAGen::openDBs ( )
private

Definition at line 225 of file CORSIKAGen_module.cc.

References allTimeWatchdog::endl, fdb, modify_metadata_with_upmu_sample_type::flist, cet::getenv(), MECModelEnuComparisons::i, ProjMan::ifdh, makeTrainCVSamples::int, LOG_DEBUG, LOG_INFO, m_fcl_ShowerInputFiles, m_IFDH, m_ShowerInputs, make_pair(), path, submit_syst::pattern, rnd_, genie::RandomGen::RndNum(), and string.

Referenced by CORSIKAGen().

225  {
226  // choose files based on m_fcl_ShowerInputFiles, copy them with ifdh, open them
227  // for c2: statement is unused
228  // sqlite3_stmt *statement;
229  // setup ifdh object
230  if (!m_IFDH) m_IFDH = new ifdh_ns::ifdh;
231  const char *ifdh_debug_env = std::getenv("IFDH_DEBUG_LEVEL");
232  if (ifdh_debug_env) {
233  LOG_INFO("CORSIKAGen") << "IFDH_DEBUG_LEVEL: " << ifdh_debug_env << "\n";
234  m_IFDH->set_debug(ifdh_debug_env);
235  }
236  // get ifdh path for each file in m_fcl_ShowerInputFiles, put into
237  // selectedflist if 1 file returned, use that file if >1 file returned,
238  // randomly select one file if 0 returned, make exeption for missing files
239  std::vector<std::pair<std::string, long>> selectedflist;
240  for (int i = 0; i < m_ShowerInputs; i++) {
241  if (m_fcl_ShowerInputFiles[i].find("*") == std::string::npos) {
242  // if there are no wildcards, don't call findMatchingFiles
243  selectedflist.push_back(std::make_pair(m_fcl_ShowerInputFiles[i], 0));
244  LOG_INFO("CorsikaGen") << "Selected" << selectedflist.back().first << "\n";
245  } else {
246  // use findMatchingFiles
247  std::vector<std::pair<std::string, long>> flist;
248  std::string path(gSystem->DirName(m_fcl_ShowerInputFiles[i].c_str()));
249  std::string pattern(gSystem->BaseName(m_fcl_ShowerInputFiles[i].c_str()));
250  flist = m_IFDH->findMatchingFiles(path, pattern);
251  unsigned int selIndex = -1;
252  if (flist.size() == 1) { // 0th element is the search path:pattern
253  selIndex = 0;
254  } else if (flist.size() > 1) {
255  double randomNumber = rnd_->RndNum().Rndm();
256  selIndex = (unsigned int)(randomNumber * (flist.size() - 1) + 0.5); // rnd with rounding, dont allow picking the 0th element
257  } else {
258  throw cet::exception("CORSIKAGen") << "No files returned for path:pattern: " << path << ":" << pattern << std::endl;
259  }
260  selectedflist.push_back(flist[selIndex]);
261  LOG_INFO("CorsikaGen")
262  << "For " << m_fcl_ShowerInputFiles[i] << ":" << pattern << "\nFound "
263  << flist.size() << " candidate files"
264  << "\nChoosing file number " << selIndex << "\n"
265  << "\nSelected " << selectedflist.back().first << "\n";
266  }
267  }
268  // do the fetching, store local filepaths in locallist
269  std::vector<std::string> locallist;
270  for (unsigned int i = 0; i < selectedflist.size(); i++) {
271  LOG_INFO("CorsikaGen") << "Fetching: " << selectedflist[i].first << " " << selectedflist[i].second << "\n";
272  std::string fetchedfile(m_IFDH->fetchInput(selectedflist[i].first));
273  LOG_DEBUG("CorsikaGen") << "Fetched; local path: " << fetchedfile;
274  locallist.push_back(fetchedfile);
275  }
276  // open the files in m_fcl_ShowerInputFilesLocalPaths with sqlite3
277  for (unsigned int i = 0; i < locallist.size(); i++) {
278  // prepare and execute statement to attach db file
279  int res = sqlite3_open(locallist[i].c_str(), &fdb[i]);
280  if (res != SQLITE_OK)
281  throw cet::exception("CORSIKAGen")
282  << "Error opening db: (" << locallist[i].c_str() << ") (" << res
283  << "): " << sqlite3_errmsg(fdb[i]) << "; memory used:<<"
284  << sqlite3_memory_used() << "/" << sqlite3_memory_highwater(0)
285  << "\n";
286  else
287  LOG_INFO("CORSIKAGen") << "Attached db " << locallist[i] << "\n";
288  }
289 
290  /* A short description of the database file:
291  Format: SQLite
292  N Table: 3
293  - input:
294  + runnr int
295  + version float
296  + nshow int
297  + model_high int
298  + model_low int
299  + eslope float
300  + erange_high float
301  + erange_low float
302  + ecuts_hadron float
303  + ecuts_muon float
304  + ecuts_electron float
305  + ecuts_photon float
306  - particles:
307  + shower int
308  + pdg int
309  + px float
310  + py float
311  + pz float
312  + x float
313  + z float
314  + t float
315  + e float
316  - showers:
317  + id int
318  + nparticles int
319 
320  */
321 }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
genie::RandomGen * rnd_
Random generator from GENIE.
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
std::vector< std::string > m_fcl_ShowerInputFiles
Set of CORSIKA shower data files to use.
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
std::string getenv(std::string const &name)
int m_ShowerInputs
Number of shower inputs to process from.
ifdh
Definition: ProjMan.py:8
const std::string path
Definition: plot_BEN.C:43
TRandom3 & RndNum(void) const
rnd number generator used by MC integrators & other numerical methods
Definition: RandomGen.h:78
ifdh_ns::ifdh * m_IFDH
(optional) flux file handling
#define LOG_INFO(stream)
Definition: Messenger.h:144
enum BeamMode string
void evgen::CORSIKAGen::populateNShowers ( )
private

Definition at line 358 of file CORSIKAGen_module.cc.

References geo::GeometryBase::DetectorBigBox(), fdb, geom(), MECModelEnuComparisons::i, makeTrainCVSamples::int, LOG_INFO, m_fcl_RandomXZShift, m_fcl_SampleTime, m_fcl_ShowerAreaExtension, m_fcl_ShowerFluxConstants, m_fcl_Toffset, m_MaxShowers, m_NShowersPerEvent, M_PI, m_ShowerBounds, m_ShowerInputs, cet::pow(), and string.

Referenced by CORSIKAGen().

358  {
359  // populate vector of the number of showers per event based on:
360  // AREA the showers are being distributed over
361  // TIME of the event (m_fcl_SampleTime)
362  // flux constants that determine the overall normalizations
363  // (m_fcl_ShowerFluxConstants) Energy range over which the sample was generated
364  // (ERANGE_*) power spectrum over which the sample was generated (ESLOPE)
365  // compute shower area based on the maximal x,z dimensions of cryostat
366  // boundaries + m_fcl_ShowerAreaExtension
368 
369  double xlo_cm = 0;
370  double xhi_cm = 0;
371  double ylo_cm = 0;
372  double yhi_cm = 0;
373  double zlo_cm = 0;
374  double zhi_cm = 0;
375 
376  geom->DetectorBigBox(&xlo_cm, &xhi_cm, &ylo_cm, &yhi_cm, &zlo_cm, &zhi_cm);
377  // add on m_fcl_ShowerAreaExtension without being clever
382  double showersArea = (m_ShowerBounds[1] / 100 - m_ShowerBounds[0] / 100) * (m_ShowerBounds[5] / 100 - m_ShowerBounds[4] / 100);
383  LOG_INFO("CORSIKAGen")
384  << "Area extended by : " << m_fcl_ShowerAreaExtension
385  << "\nShowers to be distributed betweeen: x=" << m_ShowerBounds[0] << ","
386  << m_ShowerBounds[1] << " & z=" << m_ShowerBounds[4] << ","
387  << m_ShowerBounds[5]
388  << "\nShowers to be distributed with random XZ shift: "
389  << m_fcl_RandomXZShift << " cm"
390  << "\nShowers to be distributed over area: " << showersArea << " m^2"
391  << "\nShowers to be distributed over time: " << m_fcl_SampleTime << " s"
392  << "\nShowers to be distributed with time offset: " << m_fcl_Toffset
393  << " s"
394  << "\nShowers to be distributed at y: " << m_ShowerBounds[3] << " cm";
395  // db variables
396  sqlite3_stmt *statement;
397  const std::string kStatement("select erange_high,erange_low,eslope,nshow from input");
398  double upperLimitOfEnergyRange = 0., lowerLimitOfEnergyRange = 0., energySlope = 0., oneMinusGamma = 0., EiToOneMinusGamma = 0., EfToOneMinusGamma = 0.;
399 
400  for (int i = 0; i < m_ShowerInputs; i++) {
401  // build and do query to get run info from databases
402  if (sqlite3_prepare(fdb[i], kStatement.c_str(), -1, &statement, 0) == SQLITE_OK) {
403  int res = 0;
404  res = sqlite3_step(statement);
405  if (res == SQLITE_ROW) {
406  upperLimitOfEnergyRange = sqlite3_column_double(statement, 0);
407  lowerLimitOfEnergyRange = sqlite3_column_double(statement, 1);
408  energySlope = sqlite3_column_double(statement, 2);
409  m_MaxShowers.push_back(sqlite3_column_int(statement, 3));
410  oneMinusGamma = 1 + energySlope;
411  EiToOneMinusGamma = pow(lowerLimitOfEnergyRange, oneMinusGamma);
412  EfToOneMinusGamma = pow(upperLimitOfEnergyRange, oneMinusGamma);
413  mf::LogVerbatim("CORSIKAGen")
414  << "For showers input " << i
415  << " found e_hi=" << upperLimitOfEnergyRange
416  << ", e_lo=" << lowerLimitOfEnergyRange << ", slope=" << energySlope
417  << ", k=" << m_fcl_ShowerFluxConstants[i] << "\n";
418  } else {
419  throw cet::exception("CORSIKAGen")
420  << "Unexpected sqlite3_step return value: (" << res << "); "
421  << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
422  }
423  } else {
424  throw cet::exception("CORSIKAGen")
425  << "Error preparing statement: (" << kStatement << "); "
426  << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
427  }
428  // this is computed, how?
429  double NShowers = (M_PI * showersArea * m_fcl_ShowerFluxConstants[i] * (EfToOneMinusGamma - EiToOneMinusGamma) / oneMinusGamma) * m_fcl_SampleTime;
430  m_NShowersPerEvent.push_back(NShowers);
431  mf::LogVerbatim("CORSIKAGen") << "For showers input " << i << " the number of showers per event is " << (long int)NShowers << "\n";
432  }
433 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
double m_fcl_RandomXZShift
showers won&#39;t repeatedly sample the same space [cm]
std::vector< int > m_MaxShowers
constexpr T pow(T x)
Definition: pow.h:75
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
#define M_PI
Definition: SbMath.h:34
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
std::vector< double > m_NShowersPerEvent
duration m_fcl_SampleTime; one per showerinput
double m_fcl_Toffset
Time offset of sample, defaults to zero (no offset) [s].
int m_ShowerInputs
Number of shower inputs to process from.
double m_fcl_SampleTime
Duration of sample [s].
void geom(int which=0)
Definition: geom.C:163
std::vector< double > m_fcl_ShowerFluxConstants
with each shower data file
#define LOG_INFO(stream)
Definition: Messenger.h:144
double m_fcl_ShowerAreaExtension
(e.g. 1000 will extend 10 m in -x, +x, -z, and +z) [cm]
void DetectorBigBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
enum BeamMode string
void evgen::CORSIKAGen::populateTOffset ( )
private

Definition at line 334 of file CORSIKAGen_module.cc.

References fdb, MECModelEnuComparisons::i, LOG_INFO, m_ShowerInputs, m_Toffset_corsika, string, and confusionMatrixTree::t.

Referenced by CORSIKAGen().

334  {
335  // populate TOffset_corsika by finding minimum ParticleTime from db file
336  sqlite3_stmt *statement;
337  const std::string kStatement("select min(t) from particles");
338  double t = 0.;
339  for (int i = 0; i < m_ShowerInputs; i++) {
340  // build and do query to get run min(t) from each db
341  if (sqlite3_prepare(fdb[i], kStatement.c_str(), -1, &statement, 0) == SQLITE_OK) {
342  int res = 0;
343  res = sqlite3_step(statement);
344  if (res == SQLITE_ROW) {
345  t = sqlite3_column_double(statement, 0);
346  LOG_INFO("CORSIKAGen") << "For showers input " << i << " found particles.min(t)=" << t << "\n";
347  if (i == 0 || t < m_Toffset_corsika) m_Toffset_corsika = t;
348  } else {
349  throw cet::exception("CORSIKAGen") << "Unexpected sqlite3_step return value: (" << res << "); " << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
350  }
351  } else {
352  throw cet::exception("CORSIKAGen") << "Error preparing statement: (" << kStatement << "); " << "ERROR:" << sqlite3_errmsg(fdb[i]) << "\n";
353  }
354  }
355  LOG_INFO("CORSIKAGen") << "Found corsika timeoffset [ns]: " << m_Toffset_corsika << "\n";
356 }
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
double m_Toffset_corsika
atmosphere, populated from db (optional) flux file handling
int m_ShowerInputs
Number of shower inputs to process from.
#define LOG_INFO(stream)
Definition: Messenger.h:144
enum BeamMode string
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
void evgen::CORSIKAGen::produce ( art::Event evt)
overridevirtual

Implements art::EDProducer.

Definition at line 589 of file CORSIKAGen_module.cc.

References abs(), simb::MCTruth::Add(), DetectorBigBoxCut(), simb::MCParticle::E(), energy, geom(), simb::MCTruth::GetParticle(), GetSample(), simb::kCosmicRay, LOG_INFO, m_fcl_EnhanceNNbarBkgd, m_fcl_IsBigBoxUsed, m_fcl_NNbarBkgdEnergyCut, simb::MCTruth::NParticles(), make_root_from_grid_output::pdg, simb::MCParticle::PdgCode(), art::Event::put(), and simb::MCTruth::SetOrigin().

589  {
590  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
591 
593 
594  simb::MCTruth pretruth;
595  pretruth.SetOrigin(simb::kCosmicRay);
596  GetSample(pretruth);
597  LOG_INFO("CORSIKAGen") << "GetSample() number of particles returned: " << pretruth.NParticles() << "\n";
598 
599  // DetectorBigBox cut
600  if (m_fcl_IsBigBoxUsed) this->DetectorBigBoxCut(pretruth);
601 
602  // LOG_INFO("CORSIKAGen") << "Number of particles from GetSample() crossing DetectorBigBox: " << pretruth.NParticles() << "\n";
603  // unsigned int countInTime = 0;
604  // for (unsigned int iPart = 0; iPart < pretruth.NParticles(); ++iPart) {
605  // auto particle = pretruth.GetParticle(iPart);
606  // int pdg = abs(particle.PdgCode());
607  // double px = particle.Px();
608  // double py = particle.Py();
609  // double pz = particle.Pz();
610  // double vx = particle.EndX();
611  // double vy = particle.EndY();
612  // double vz = particle.EndZ();
613  // double vt = particle.EndT();
614 
615  // std::cout << Form("PRETRUTH (%05i): ", iPart) << std::fixed
616  // << std::setw(10) << pdg << ", " << std::setw(15)
617  // << std::setprecision(5) << vx << ", " << std::setw(15)
618  // << std::setprecision(5) << vy << ", " << std::setw(15)
619  // << std::setprecision(5) << vz << ", " << std::setw(15)
620  // << std::setprecision(5) << vt / 1000. << " us, "
621  // << std::setw(15) << std::setprecision(5) << px << ", "
622  // << std::setw(15) << std::setprecision(5) << py << ", "
623  // << std::setw(15) << std::setprecision(5) << pz << std::endl;
624 
625  // if (vt / 1000. <= 0 || vt / 1000. >= 550) {
626  // continue;
627  // } else {
628  // countInTime++;
629  // }
630  // }
631  // LOG_INFO("CORSIKAGen") << "Number of particles from GetSample() crossing DetectorBigBox and in time: " << countInTime << "\n";
632 
634  simb::MCTruth pretruth_nnbar_bkgd;
635  for (unsigned int iPart = 0; iPart < pretruth.NParticles(); ++iPart) {
636  simb::MCParticle particle = pretruth.GetParticle(iPart);
637  int pdg = abs(particle.PdgCode());
638  double energy = particle.E();
639  if (pdg == 2112 || pdg == -2112 || pdg == 22) {
640  if (energy >= m_fcl_NNbarBkgdEnergyCut) pretruth_nnbar_bkgd.Add(particle);
641  }
642  }
643  truthcol->push_back(pretruth_nnbar_bkgd);
644  } else {
645  truthcol->push_back(pretruth);
646  }
647  evt.put(std::move(truthcol));
648 
649  return;
650 }
double E(const int i=0) const
Definition: MCParticle.h:232
bool m_fcl_EnhanceNNbarBkgd
Generating only cosmogenic (anti)neutrons and gammas if TRUE.
int PdgCode() const
Definition: MCParticle.h:211
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
void abs(TH1 *hist)
int NParticles() const
Definition: MCTruth.h:74
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void DetectorBigBoxCut(simb::MCTruth &truth)
double energy
Definition: plottest35.C:25
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
double m_fcl_NNbarBkgdEnergyCut
Energy (in GeV) cut for NNbar background particles (neutron, antineutron, gamma)
void geom(int which=0)
Definition: geom.C:163
Event generator information.
Definition: MCTruth.h:32
#define LOG_INFO(stream)
Definition: Messenger.h:144
void GetSample(simb::MCTruth &)
Cosmic rays.
Definition: MCTruth.h:24
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

Definition at line 155 of file CORSIKAGen_module.cc.

References d, dx, dy, dz, and MECModelEnuComparisons::i.

Referenced by GetSample().

160  {
161  // we want to project backwards, so take mirror of momentum
162  const double dxyz[3] = {-indxyz[0], -indxyz[1], -indxyz[2]};
163  // Compute the distances to the x/y/z walls
164  double dx = 99.E99;
165  double dy = 99.E99;
166  double dz = 99.E99;
167  if (dxyz[0] > 0.0) {
168  dx = (xhi - xyz[0]) / dxyz[0];
169  } else if (dxyz[0] < 0.0) {
170  dx = (xlo - xyz[0]) / dxyz[0];
171  }
172  if (dxyz[1] > 0.0) {
173  dy = (yhi - xyz[1]) / dxyz[1];
174  } else if (dxyz[1] < 0.0) {
175  dy = (ylo - xyz[1]) / dxyz[1];
176  }
177  if (dxyz[2] > 0.0) {
178  dz = (zhi - xyz[2]) / dxyz[2];
179  } else if (dxyz[2] < 0.0) {
180  dz = (zlo - xyz[2]) / dxyz[2];
181  }
182  // Choose the shortest distance
183  double d = 0.0;
184  if (dx < dy && dx < dz)
185  d = dx;
186  else if (dy < dz && dy < dx)
187  d = dy;
188  else if (dz < dx && dz < dy)
189  d = dz;
190  // Make the step
191  for (int i = 0; i < 3; ++i) {
192  xyzout[i] = xyz[i] + dxyz[i] * d;
193  }
194 }
double dy[NP][NC]
double dx[NP][NC]
double dz[NP][NC]
Float_t d
Definition: plot.C:236
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Referenced by art::RootOutput::endJob().

void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited
double evgen::CORSIKAGen::wrapvar ( const double  var,
const double  low,
const double  high 
)
private

Definition at line 323 of file CORSIKAGen_module.cc.

References stan::math::floor().

Referenced by GetSample().

323  {
324  // wrap variable so that it's always between low and high
325  return (var - (high - low) * floor(var / (high - low))) + low;
326 }
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
double evgen::CORSIKAGen::wrapvarBoxNo ( const double  var,
const double  low,
const double  high,
int boxno 
)
private

Definition at line 328 of file CORSIKAGen_module.cc.

References stan::math::floor(), and makeTrainCVSamples::int.

Referenced by GetSample().

328  {
329  // wrap variable so that it's always between low and high
330  boxno = int(floor(var / (high - low)));
331  return (var - (high - low) * floor(var / (high - low))) + low;
332 }
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11

Member Data Documentation

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

Pointers to sqlite3 database object, max of 5.

Definition at line 98 of file CORSIKAGen_module.cc.

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

int evgen::CORSIKAGen::m_fcl_Cycle
private

MC cycle generation number.

Definition at line 118 of file CORSIKAGen_module.cc.

Referenced by endSubRun().

bool evgen::CORSIKAGen::m_fcl_EnhanceNNbarBkgd
private

Generating only cosmogenic (anti)neutrons and gammas if TRUE.

Definition at line 115 of file CORSIKAGen_module.cc.

Referenced by produce().

bool evgen::CORSIKAGen::m_fcl_IsBigBoxUsed
private

Do we need to use the BigBox cut? The cosmic rays must go through the DetectorBigBox. Otherwise, don't store them

Definition at line 111 of file CORSIKAGen_module.cc.

Referenced by produce().

double evgen::CORSIKAGen::m_fcl_NNbarBkgdEnergyCut
private

Energy (in GeV) cut for NNbar background particles (neutron, antineutron, gamma)

Definition at line 116 of file CORSIKAGen_module.cc.

Referenced by produce().

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

Height to which particles will be projected [cm].

Definition at line 101 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), and GetSample().

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

showers won't repeatedly sample the same space [cm]

Each shower will be shifted by a random amount in xz so that

Definition at line 109 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

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

Duration of sample [s].

Definition at line 105 of file CORSIKAGen_module.cc.

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

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

(e.g. 1000 will extend 10 m in -x, +x, -z, and +z) [cm]

Extend distribution of corsika particles in x,z by this much

Definition at line 107 of file CORSIKAGen_module.cc.

Referenced by populateNShowers().

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

with each shower data file

Set of flux constants to be associated

Definition at line 103 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), and populateNShowers().

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

Set of CORSIKA shower data files to use.

Definition at line 102 of file CORSIKAGen_module.cc.

Referenced by CORSIKAGen(), and openDBs().

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

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

Definition at line 106 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

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

(optional) flux file handling

Definition at line 97 of file CORSIKAGen_module.cc.

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

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

Definition at line 90 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

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

duration m_fcl_SampleTime; one per showerinput

Number of showers to put in each event of

Definition at line 88 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

double evgen::CORSIKAGen::m_ShowerBounds[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 91 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateNShowers().

int evgen::CORSIKAGen::m_ShowerInputs = 0
private

Number of shower inputs to process from.

Definition at line 87 of file CORSIKAGen_module.cc.

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

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

atmosphere, populated from db (optional) flux file handling

Timing offset to account for propagation time through

Definition at line 95 of file CORSIKAGen_module.cc.

Referenced by GetSample(), and populateTOffset().

genie::RandomGen* evgen::CORSIKAGen::rnd_
private

Random generator from GENIE.

Definition at line 85 of file CORSIKAGen_module.cc.

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

TStopwatch evgen::CORSIKAGen::stopwatch_
private

keep track of how long it takes to run the job

Definition at line 84 of file CORSIKAGen_module.cc.

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


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