23 #include "TDatabasePDG.h" 28 #include "GENIE/Framework/Numerical/RandomGen.h" 34 #include "TStopwatch.h" 35 #include "Utilities/AssociationUtil.h" 69 double wrapvar(
const double var,
const double low,
const double high);
70 double wrapvarBoxNo(
const double var,
const double low,
const double high,
76 double xhi,
double ylo,
double yhi,
double zlo,
80 const double xlo,
const double xhi,
const double ylo,
81 const double yhi,
const double zlo,
const double zhi,
139 throw cet::exception(
"CORSIKAGen") <<
"ShowerInputFiles and ShowerFluxConstants have different or invalid sizes!" <<
"\n";
150 produces<std::vector<simb::MCTruth>>();
151 produces<sumdata::SubRunData, art::InSubRun>();
152 produces<sumdata::RunData, art::InRun>();
156 const double indxyz[],
157 const double xlo,
const double xhi,
158 const double ylo,
const double yhi,
159 const double zlo,
const double zhi,
162 const double dxyz[3] = {-indxyz[0], -indxyz[1], -indxyz[2]};
168 dx = (xhi - xyz[0]) / dxyz[0];
169 }
else if (dxyz[0] < 0.0) {
170 dx = (xlo - xyz[0]) / dxyz[0];
173 dy = (yhi - xyz[1]) / dxyz[1];
174 }
else if (dxyz[1] < 0.0) {
175 dy = (ylo - xyz[1]) / dxyz[1];
178 dz = (zhi - xyz[2]) / dxyz[2];
179 }
else if (dxyz[2] < 0.0) {
180 dz = (zlo - xyz[2]) / dxyz[2];
184 if (dx < dy && dx < dz)
186 else if (dy < dz && dy < dx)
188 else if (dz < dx && dz < dy)
191 for (
int i = 0;
i < 3; ++
i) {
192 xyzout[
i] = xyz[
i] + dxyz[
i] *
d;
198 sqlite3_close(
fdb[
i]);
205 <<
"real time to produce file: " <<
stopwatch_.RealTime();
213 run.
put(std::move(runcol));
222 subrun.
put(std::move(sd));
231 const char *ifdh_debug_env =
std::getenv(
"IFDH_DEBUG_LEVEL");
232 if (ifdh_debug_env) {
233 MF_LOG_INFO(
"CORSIKAGen") <<
"IFDH_DEBUG_LEVEL: " << ifdh_debug_env <<
"\n";
234 m_IFDH->set_debug(ifdh_debug_env);
239 std::vector<std::pair<std::string, long>> selectedflist;
244 MF_LOG_INFO(
"CorsikaGen") <<
"Selected" << selectedflist.back().first <<
"\n";
247 std::vector<std::pair<std::string, long>>
flist;
250 flist =
m_IFDH->findMatchingFiles(path, pattern);
251 unsigned int selIndex = -1;
252 if (flist.size() == 1) {
254 }
else if (flist.size() > 1) {
256 selIndex = (
unsigned int)(randomNumber * (flist.size() - 1) + 0.5);
258 throw cet::exception(
"CORSIKAGen") <<
"No files returned for path:pattern: " << path <<
":" << pattern <<
std::endl;
260 selectedflist.push_back(flist[selIndex]);
263 << flist.size() <<
" candidate files" 264 <<
"\nChoosing file number " << selIndex <<
"\n" 265 <<
"\nSelected " << selectedflist.back().first <<
"\n";
269 std::vector<std::string> locallist;
270 for (
unsigned int i = 0;
i < selectedflist.size();
i++) {
271 MF_LOG_INFO(
"CorsikaGen") <<
"Fetching: " << selectedflist[
i].first <<
" " << selectedflist[
i].second <<
"\n";
273 MF_LOG_DEBUG(
"CorsikaGen") <<
"Fetched; local path: " << fetchedfile;
274 locallist.push_back(fetchedfile);
277 for (
unsigned int i = 0;
i < locallist.size();
i++) {
279 int res = sqlite3_open(locallist[
i].c_str(), &
fdb[
i]);
280 if (res != SQLITE_OK)
282 <<
"Error opening db: (" << locallist[
i].c_str() <<
") (" << res
283 <<
"): " << sqlite3_errmsg(
fdb[i]) <<
"; memory used:<<" 284 << sqlite3_memory_used() <<
"/" << sqlite3_memory_highwater(0)
287 MF_LOG_INFO(
"CORSIKAGen") <<
"Attached db " << locallist[
i] <<
"\n";
325 return (var - (high - low) * floor(var / (high - low))) + low;
330 boxno =
int(floor(var / (high - low)));
331 return (var - (high - low) * floor(var / (high - low))) + low;
336 sqlite3_stmt *statement;
337 const std::string kStatement(
"select min(t) from particles");
341 if (sqlite3_prepare(
fdb[
i], kStatement.c_str(), -1, &statement, 0) == SQLITE_OK) {
343 res = sqlite3_step(statement);
344 if (res == SQLITE_ROW) {
345 t = sqlite3_column_double(statement, 0);
346 MF_LOG_INFO(
"CORSIKAGen") <<
"For showers input " << i <<
" found particles.min(t)=" << t <<
"\n";
349 throw cet::exception(
"CORSIKAGen") <<
"Unexpected sqlite3_step return value: (" << res <<
"); " <<
"ERROR:" << sqlite3_errmsg(
fdb[i]) <<
"\n";
352 throw cet::exception(
"CORSIKAGen") <<
"Error preparing statement: (" << kStatement <<
"); " <<
"ERROR:" << sqlite3_errmsg(
fdb[i]) <<
"\n";
376 geom->
DetectorBigBox(&xlo_cm, &xhi_cm, &ylo_cm, &yhi_cm, &zlo_cm, &zhi_cm);
384 <<
"Area extended by : " << m_fcl_ShowerAreaExtension
385 <<
"\nShowers to be distributed betweeen: x=" <<
m_ShowerBounds[0] <<
"," 388 <<
"\nShowers to be distributed with random XZ shift: " 390 <<
"\nShowers to be distributed over area: " << showersArea <<
" m^2" 392 <<
"\nShowers to be distributed with time offset: " <<
m_fcl_Toffset 394 <<
"\nShowers to be distributed at y: " <<
m_ShowerBounds[3] <<
" cm";
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.;
402 if (sqlite3_prepare(
fdb[
i], kStatement.c_str(), -1, &statement, 0) == SQLITE_OK) {
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);
414 <<
"For showers input " << i
415 <<
" found e_hi=" << upperLimitOfEnergyRange
416 <<
", e_lo=" << lowerLimitOfEnergyRange <<
", slope=" << energySlope
420 <<
"Unexpected sqlite3_step return value: (" << res <<
"); " 421 <<
"ERROR:" << sqlite3_errmsg(
fdb[i]) <<
"\n";
425 <<
"Error preparing statement: (" << kStatement <<
"); " 426 <<
"ERROR:" << sqlite3_errmsg(
fdb[i]) <<
"\n";
431 mf::LogVerbatim(
"CORSIKAGen") <<
"For showers input " << i <<
" the number of showers per event is " << (
long int)NShowers <<
"\n";
451 static TDatabasePDG *pdgt = TDatabasePDG::Instance();
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)");
467 geom->
WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);
470 double fBoxDelta = 1.e-5;
482 long int nShowerCntr = 0;
485 double px, py, pz,
x,
z, tParticleTime, etot, showerTime = 0., showerTimex = 0., showerTimez = 0., showerXOffset = 0., showerZOffset = 0.,
t;
489 <<
". Maximum number of showers in database: " <<
m_MaxShowers[
i];
490 while (nShowerCntr > 0) {
492 if (nShowerCntr > m_MaxShowers[i]) nShowerQry = m_MaxShowers[
i];
493 else nShowerQry = nShowerCntr;
497 TString kthisStatement =
TString::Format(kStatement.Data(), thisrnd, nShowerQry, thisrnd);
498 MF_LOG_INFO(
"CORSIKAGen") <<
"Executing: " << kthisStatement;
499 if (sqlite3_prepare(
fdb[i], kthisStatement.Data(), -1, &statement, 0) == SQLITE_OK) {
503 res = sqlite3_step(statement);
504 if (res == SQLITE_ROW) {
517 shower = sqlite3_column_int(statement, 0);
518 if (shower != lastShower) {
528 pdg = sqlite3_column_int(statement, 1);
531 TParticlePDG *pdgp = pdgt->GetParticle(pdg);
532 if (pdgp) m = pdgp->Mass();
536 px = sqlite3_column_double(statement, 4);
537 py = sqlite3_column_double(statement, 3);
538 pz = -sqlite3_column_double(statement, 2);
539 etot = sqlite3_column_double(statement, 8);
542 int boxnoX = 0, boxnoZ = 0;
545 tParticleTime = sqlite3_column_double(statement, 7);
567 double dx[3] = {px, py, pz};
569 TLorentzVector
pos(xyzo[0], xyzo[1], xyzo[2], t);
570 TLorentzVector mom(px, py, pz, etot);
575 }
else if (res == SQLITE_DONE) {
578 throw cet::exception(
"CORSIKAGen") <<
"Unexpected sqlite3_step return value: (" << res <<
"); " <<
"ERROR:" << sqlite3_errmsg(
fdb[i]) <<
"\n";
582 throw cet::exception(
"CORSIKAGen") <<
"Error preparing statement: (" << kthisStatement <<
"); " <<
"ERROR:" << sqlite3_errmsg(
fdb[i]) <<
"\n";
584 nShowerCntr = nShowerCntr - nShowerQry;
590 std::unique_ptr<std::vector<simb::MCTruth>> truthcol(
new std::vector<simb::MCTruth>);
597 MF_LOG_INFO(
"CORSIKAGen") <<
"GetSample() number of particles returned: " << pretruth.
NParticles() <<
"\n";
635 for (
unsigned int iPart = 0; iPart < pretruth.
NParticles(); ++iPart) {
639 if (pdg == 2112 || pdg == -2112 || pdg == 22) {
643 truthcol->push_back(pretruth_nnbar_bkgd);
645 truthcol->push_back(pretruth);
647 evt.
put(std::move(truthcol));
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";
684 double px = particle.
Px() / particle.
P();
685 double py = particle.
Py() / particle.
P();
686 double pz = particle.
Pz() / particle.
P();
688 double vx = particle.
EndX();
689 double vy = particle.
EndY();
690 double vz = particle.
EndZ();
702 double xyz[3] = {vx, vy, vz};
703 double dxyz[3] = {px, py, pz};
709 part.AddTrajectoryPoint(TLorentzVector(vx, vy, vz, particle.
T()), particle.
Momentum());
719 const double dxyz[],
double xlo,
721 double zlo,
double zhi) {
723 if (xyz[0] >= xlo && xyz[0] <= xhi && xyz[1] >= ylo && xyz[1] <= yhi && xyz[2] >= zlo && xyz[2] <= zhi)
727 double dx = 0.,
dy = 0.,
dz = 0.;
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;
740 if (x >= xlo && x <= xhi && z >= zlo && z <= zhi)
return true;
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;
753 if (x >= xlo && x <= xhi && z >= zlo && z <= zhi)
return true;
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;
766 if (x >= xlo && x <= xhi && y >= ylo && y <= yhi)
return true;
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;
779 if (x >= xlo && x <= xhi && y >= ylo && y <= yhi)
return true;
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;
792 if (z >= zlo && z <= zhi && y >= ylo && y <= yhi)
return true;
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;
805 if (z >= zlo && z <= zhi && y >= ylo && y <= yhi)
return true;
double E(const int i=0) const
bool m_fcl_EnhanceNNbarBkgd
Generating only cosmogenic (anti)neutrons and gammas if TRUE.
double wrapvar(const double var, const double low, const double high)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
void produce(art::Event &evt) override
double Py(const int i=0) const
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
double m_fcl_ProjectToHeight
Height to which particles will be projected [cm].
Float_t y1[n_points_granero]
void SetOrigin(simb::Origin_t origin)
static RandomGen * Instance()
Access instance.
Float_t x1[n_points_granero]
double m_fcl_RandomXZShift
showers won't repeatedly sample the same space [cm]
std::vector< int > m_MaxShowers
double Px(const int i=0) const
genie::RandomGen * rnd_
Random generator from GENIE.
::xsd::cxx::tree::exception< char > exception
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
DEFINE_ART_MODULE(TestTMapFile)
std::string Process() const
std::vector< std::string > m_fcl_ShowerInputFiles
Set of CORSIKA shower data files to use.
::xsd::cxx::tree::time< char, simple_type > time
A singleton holding random number generator classes. All random number generation in GENIE should tak...
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
virtual void endSubRun(art::SubRun &run)
void WorldBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
std::string getenv(std::string const &name)
TStopwatch stopwatch_
keep track of how long it takes to run the job
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 P(const int i=0) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
double m_Toffset_corsika
atmosphere, populated from db (optional) flux file handling
void DetectorBigBoxCut(simb::MCTruth &truth)
double T(const int i=0) const
Near Detector in the NuMI cavern.
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)
#define MF_LOG_INFO(category)
const simb::MCParticle & GetParticle(int i) const
TRandom3 & RndNum(void) const
rnd number generator used by MC integrators & other numerical methods
CORSIKAGen(fhicl::ParameterSet const &p)
ifdh_ns::ifdh * m_IFDH
(optional) flux file handling
void Add(simb::MCParticle const &part)
double m_fcl_SampleTime
Duration of sample [s].
double m_fcl_NNbarBkgdEnergyCut
Energy (in GeV) cut for NNbar background particles (neutron, antineutron, gamma)
const TLorentzVector & Momentum(const int i=0) const
double Pz(const int i=0) const
std::vector< double > m_fcl_ShowerFluxConstants
with each shower data file
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
virtual void beginRun(art::Run &run)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
bool isIntersectTheBox(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
virtual void beginSubRun(art::SubRun &run)
Event generator information.
double m_fcl_ShowerAreaExtension
(e.g. 1000 will extend 10 m in -x, +x, -z, and +z) [cm]
void GetSample(simb::MCTruth &)
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Module to generate only pions from cosmic rays.
Encapsulate the geometry of one entire detector (near, far, ndos)
int m_fcl_Cycle
MC cycle generation number.
void DetectorBigBox(double *xlo_cm, double *xhi_cm, double *ylo_cm, double *yhi_cm, double *zlo_cm, double *zhi_cm) const
std::string FileBaseName() const
void SetSeed(long int seed)