CORSIKAGen_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CORSIKAGen
3 // Plugin Type: producer (art v2_13_00)
4 // File: CORSIKAGen_module.cc
5 //
6 // Generated at Thu Jan 30 15:45:47 2020 by Dung Phan using cetskelgen
7 // from cetlib version v3_06_01.
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include <sqlite3.h>
11 #include <unistd.h>
12 
13 #include <cassert>
14 #include <cstdlib>
15 #include <ctime>
16 #include <map>
17 #include <memory>
18 #include <sstream>
19 #include <string>
20 #include <iomanip>
21 #include <vector>
22 
23 #include "TDatabasePDG.h"
24 #include "TString.h"
25 #include "TSystem.h"
26 
27 // GENIE includes
28 #include "GENIE/Framework/Numerical/RandomGen.h"
29 #include "Geometry/Geometry.h"
30 #include "SummaryData/POTSum.h"
31 #include "SummaryData/RunData.h"
32 #include "SummaryData/SpillData.h"
33 #include "SummaryData/SubRunData.h"
34 #include "TStopwatch.h"
35 #include "Utilities/AssociationUtil.h"
44 #include "fhiclcpp/ParameterSet.h"
45 #include "ifdh.h" //to handle flux files
49 
50 namespace evgen {
51 class CORSIKAGen;
52 }
53 
55  public:
56  explicit CORSIKAGen(fhicl::ParameterSet const &p);
57  virtual ~CORSIKAGen();
58 
59  void produce(art::Event &evt) override;
60  virtual void beginRun(art::Run &run);
61  virtual void beginSubRun(art::SubRun& run);
62  virtual void endSubRun(art::SubRun& run);
63 
64  private:
65  void openDBs();
66  void populateNShowers();
67  void populateTOffset();
68  void GetSample(simb::MCTruth &);
69  double wrapvar(const double var, const double low, const double high);
70  double wrapvarBoxNo(const double var, const double low, const double high,
71  int &boxno);
72 
73  void DetectorBigBoxCut(simb::MCTruth &truth);
74 
75  bool isIntersectTheBox(const double xyz[], const double dxyz[], double xlo,
76  double xhi, double ylo, double yhi, double zlo,
77  double zhi);
78 
79  void ProjectToBoxEdge(const double xyz[], const double dxyz[],
80  const double xlo, const double xhi, const double ylo,
81  const double yhi, const double zlo, const double zhi,
82  double xyzout[]);
83 
84  TStopwatch stopwatch_; ///< keep track of how long it takes to run the job
85  genie::RandomGen *rnd_; ///< Random generator from GENIE
86 
87  int m_ShowerInputs = 0; ///< Number of shower inputs to process from
88  std::vector<double> m_NShowersPerEvent; ///< Number of showers to put in each event of
89  ///< duration m_fcl_SampleTime; one per showerinput
90  std::vector<int> m_MaxShowers; //< Max number of showers to query, one per showerinput
91  double m_ShowerBounds[6] = { 0., 0.,
92  0., 0.,
93  0., 0.}; ///< Boundaries of area over which showers are to be distributed
94  ///< (x(min), x(max), _unused_, y, z(min), z(max) )
95  double m_Toffset_corsika = 0.; ///< Timing offset to account for propagation time through
96  ///< atmosphere, populated from db (optional) flux file handling
97  ifdh_ns::ifdh *m_IFDH = 0; ///< (optional) flux file handling
98  sqlite3 *fdb[5]; ///< Pointers to sqlite3 database object, max of 5
99 
100  // fcl parameters
101  double m_fcl_ProjectToHeight = 0.; ///< Height to which particles will be projected [cm]
102  std::vector<std::string> m_fcl_ShowerInputFiles; ///< Set of CORSIKA shower data files to use
103  std::vector<double> m_fcl_ShowerFluxConstants; ///< Set of flux constants to be associated
104  ///< with each shower data file
105  double m_fcl_SampleTime = 0.; ///< Duration of sample [s]
106  double m_fcl_Toffset = 0.; ///< Time offset of sample, defaults to zero (no offset) [s]
107  double m_fcl_ShowerAreaExtension = 0.; ///< Extend distribution of corsika particles in x,z by this much
108  ///< (e.g. 1000 will extend 10 m in -x, +x, -z, and +z) [cm]
109  double m_fcl_RandomXZShift = 0.; ///< Each shower will be shifted by a random amount in xz so that
110  ///< showers won't repeatedly sample the same space [cm]
111  bool m_fcl_IsBigBoxUsed; ///< Do we need to use the BigBox cut? The cosmic
112  ///< rays must go through the DetectorBigBox.
113  ///< Otherwise, don't store them
114 
115  bool m_fcl_EnhanceNNbarBkgd; ///< Generating only cosmogenic (anti)neutrons and gammas if TRUE
116  double m_fcl_NNbarBkgdEnergyCut; ///< Energy (in GeV) cut for NNbar background particles (neutron, antineutron, gamma)
117 
118  int m_fcl_Cycle; ///< MC cycle generation number.
119 };
120 
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 }
154 
155 void evgen::CORSIKAGen::ProjectToBoxEdge(const double xyz[],
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,
160  double xyzout[]) {
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 }
195 
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 }
207 
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 }
215 
217 }
218 
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 }
224 
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 }
322 
323 double evgen::CORSIKAGen::wrapvar(const double var, const double low, const double high) {
324  // wrap variable so that it's always between low and high
325  return (var - (high - low) * floor(var / (high - low))) + low;
326 }
327 
328 double evgen::CORSIKAGen::wrapvarBoxNo(const double var, const double low, const double high, int &boxno) {
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 }
333 
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 }
357 
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 }
434 
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 }
588 
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 }
651 
652 ////////////////////////////////////////////////////////////////////////
653 // Method to cut the events not going through the DetectorBigBox,
654 // which is defined as the box surrounding the DetectorEnclosureBox
655 // by adding an additional length in each of the dimensions.
656 // The additional length is defined in detectorbigbox.xml
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 }
717 
718 bool evgen::CORSIKAGen::isIntersectTheBox(const double xyz[],
719  const double dxyz[], double xlo,
720  double xhi, double ylo, double yhi,
721  double zlo, double zhi) {
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 }
810 
double E(const int i=0) const
Definition: MCParticle.h:232
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
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
int PdgCode() const
Definition: MCParticle.h:211
void produce(art::Event &evt) override
double Py(const int i=0) const
Definition: MCParticle.h:230
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:256
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double EndZ() const
Definition: MCParticle.h:227
Int_t ipart
Definition: Style.C:10
double m_fcl_ProjectToHeight
Height to which particles will be projected [cm].
Float_t y1[n_points_granero]
Definition: compare.C:5
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
Float_t x1[n_points_granero]
Definition: compare.C:5
int Mother() const
Definition: MCParticle.h:212
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
double Mass() const
Definition: MCParticle.h:238
constexpr T pow(T x)
Definition: pow.h:75
double Px(const int i=0) const
Definition: MCParticle.h:229
genie::RandomGen * rnd_
Random generator from GENIE.
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void abs(TH1 *hist)
int StatusCode() const
Definition: MCParticle.h:210
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
int NParticles() const
Definition: MCTruth.h:74
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
std::string Process() const
Definition: MCParticle.h:214
Particle class.
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
double EndY() const
Definition: MCParticle.h:226
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
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
#define M_PI
Definition: SbMath.h:34
Definition: Run.h:31
int TrackId() const
Definition: MCParticle.h:209
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double sd(Eigen::VectorXd x)
TString part[npart]
Definition: Style.C:32
Int_t npart
Definition: Style.C:7
double dy[NP][NC]
double dx[NP][NC]
sqlite3 * fdb[5]
Pointers to sqlite3 database object, max of 5.
double dz[NP][NC]
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
Definition: MCParticle.h:233
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
double m_Toffset_corsika
atmosphere, populated from db (optional) flux file handling
void DetectorBigBoxCut(simb::MCTruth &truth)
double energy
Definition: plottest35.C:25
double T(const int i=0) const
Definition: MCParticle.h:223
Near Detector in the NuMI cavern.
double m_fcl_Toffset
Time offset of sample, defaults to zero (no offset) [s].
Float_t d
Definition: plot.C:236
int m_ShowerInputs
Number of shower inputs to process from.
double wrapvarBoxNo(const double var, const double low, const double high, int &boxno)
ifdh
Definition: ProjMan.py:8
z
Definition: test.py:28
Definition: run.py:1
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
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
CORSIKAGen(fhicl::ParameterSet const &p)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
ifdh_ns::ifdh * m_IFDH
(optional) flux file handling
double m_fcl_SampleTime
Duration of sample [s].
ProductID put(std::unique_ptr< PROD > &&)
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
double m_fcl_NNbarBkgdEnergyCut
Energy (in GeV) cut for NNbar background particles (neutron, antineutron, gamma)
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
void geom(int which=0)
Definition: geom.C:163
double Pz(const int i=0) const
Definition: MCParticle.h:231
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)
Definition: Style.cxx:154
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.
Definition: MCTruth.h:32
#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]
double EndX() const
Definition: MCParticle.h:225
void GetSample(simb::MCTruth &)
Helper for AttenCurve.
Definition: Path.h:10
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.
Cosmic rays.
Definition: MCTruth.h:24
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)
Definition: RandomGen.cxx:90