LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
SingleGen_module.cc
Go to the documentation of this file.
1 
10 // C++ includes.
11 #include <cctype> // std::tolower()
12 #include <cmath>
13 #include <initializer_list>
14 #include <iterator>
15 #include <map>
16 #include <memory>
17 #include <sstream>
18 #include <string>
19 
20 // Framework includes
26 #include "cetlib/filesystem.h"
27 #include "cetlib/search_path.h"
28 #include "cetlib_except/exception.h"
29 #include "fhiclcpp/types/Atom.h"
30 #include "fhiclcpp/types/Comment.h"
31 #include "fhiclcpp/types/Name.h"
35 
36 // nurandom includes
38 
39 // nusimdata includes
42 
43 // LArSoft libraries
48 
49 #include "TDatabasePDG.h"
50 #include "TFile.h"
51 #include "TH1.h"
52 #include "TH2.h"
53 #include "TVector3.h"
54 
55 #include "CLHEP/Random/RandFlat.h"
56 #include "CLHEP/Random/RandGaussQ.h"
57 
58 namespace evgen {
59 
61  class SingleGen : public art::EDProducer {
62 
63  public:
64  struct Config {
65  using Name = fhicl::Name;
67 
69  Name("ParticleSelectionMode"),
70  Comment("generate one particle, or one particle per PDG ID: " +
72 
74  Name("PadOutVectors"),
75  Comment("if true, all per-PDG-ID quantities must contain only one value, which is then "
76  "used for all PDG IDs")};
77 
79  Comment("PDG ID of the particles to be generated; this is the key "
80  "for the other options marked as \"per PDG ID\"")};
81 
83  Name("PDist"),
84  Comment("momentum distribution type: " + presentOptions(DistributionNames)),
86 
88  Comment("central momentum (GeV/c) to generate"),
89  [this]() { return !fromHistogram(PDist()); }};
90 
92  Comment("variation in momenta (GeV/c)"),
93  [this]() { return !fromHistogram(PDist()); }};
94 
96  Name("X0"),
97  Comment("central x position (cm) in world coordinates [per PDG ID]")};
98 
100  Name("Y0"),
101  Comment("central y position (cm) in world coordinates [per PDG ID]")};
102 
104  Name("Z0"),
105  Comment("central z position (cm) in world coordinates [per PDG ID]")};
106 
107  fhicl::Sequence<double> T0{Name("T0"), Comment("central time (ns) [per PDG ID]")};
108 
110  Name("SigmaX"),
111  Comment("variation (radius or RMS) in x position (cm) [per PDG ID]")};
112 
114  Name("SigmaY"),
115  Comment("variation (radius or RMS) in y position (cm) [per PDG ID]")};
116 
118  Name("SigmaZ"),
119  Comment("variation (radius or RMS) in z position (cm) [per PDG ID]")};
120 
122  Name("SigmaT"),
123  Comment("variation (semi-interval or RMS) in time (ns) [per PDG ID]")};
124 
126  Comment("distribution of starting position: " +
128 
130  Name("TDist"),
131  Comment("time distribution type: " + presentOptions(DistributionNames, true, {kHIST}))};
132 
134  Name("SingleVertex"),
135  Comment("if true, all particles are produced at the same location"),
136  false};
137 
139  Comment("angle from Z axis on world X-Z plane (degrees)")};
140 
142  Comment("angle from Z axis on world Y-Z plane (degrees)")};
143 
145  Comment("variation in angle in X-Z plane (degrees)")};
146 
148  Comment("variation in angle in Y-Z plane (degrees)")};
149 
151  Name("AngleDist"),
152  Comment("angular distribution type: " + presentOptions(DistributionNames)),
154 
156  Name("HistogramFile"),
157  Comment("ROOT file containing the required distributions for the generation"),
158  [this]() { return fromHistogram(AngleDist()) || fromHistogram(PDist()); }};
159 
161  Name("PHist"),
162  Comment("name of the histograms of momentum distributions"),
163  [this]() { return fromHistogram(PDist()); }};
164 
165  /*
166  fhicl::Sequence<std::string> ThetaPhiHist{
167  Name("ThetaPhiHist"),
168  Comment("name of the histograms of angular (theta/phi) distribution"),
169  [this](){ return fromHistogram(AngleDist()); }
170  };
171  */
173  Name("ThetaXzYzHist"),
174  Comment("name of the histograms of angular (X-Z and Y-Z) distribution"),
175  [this]() { return fromHistogram(AngleDist()); }};
176 
178  Name("Seed"),
179  Comment("override the random number generator seed"),
181 
182  private:
184  bool fromHistogram(std::string const& key) const;
185 
186  }; // struct Config
187 
189 
190  explicit SingleGen(Parameters const& config);
191 
192  // This is called for each event.
193  void produce(art::Event& evt) override;
194 
195  private:
197  static const std::map<int, std::string> ParticleSelectionModeNames;
199  static const std::map<int, std::string> DistributionNames;
200 
201  void SampleOne(unsigned int i, simb::MCTruth& mct);
202  void SampleMany(simb::MCTruth& mct);
203  void Sample(simb::MCTruth& mct);
204  void printVecs(std::vector<std::string> const& list);
205  bool PadVector(std::vector<double>& vec);
206  double SelectFromHist(const TH1& h);
207  void SelectFromHist(const TH2& h, double& x, double& y);
209  double particle_mass(int particle_id);
211  simb::MCParticle mc_particle(int track_id, int particle_id, double mass);
212 
215 
216  static constexpr int kSelectAllParts = 0;
217  static constexpr int kSelectOneRandPart =
218  1;
219 
223 
224  static constexpr int kUNIF = 0;
225  static constexpr int kGAUS = 1;
226  static constexpr int kHIST = 2;
227 
229  int fMode;
230  bool fPadOutVectors;
233  std::vector<int> fPDG;
237  std::vector<double> fP0;
238  std::vector<double> fSigmaP;
239  int fPDist;
240  std::vector<double> fX0;
241  std::vector<double> fY0;
242  std::vector<double> fZ0;
243  std::vector<double> fT0;
244  std::vector<double> fSigmaX;
245  std::vector<double> fSigmaY;
246  std::vector<double> fSigmaZ;
247  std::vector<double> fSigmaT;
248  int fPosDist;
249  int fTDist;
251  std::vector<double> fTheta0XZ;
252  std::vector<double> fTheta0YZ;
253  std::vector<double> fSigmaThetaXZ;
254  std::vector<double> fSigmaThetaYZ;
256  std::string fHistFileName;
257  std::vector<std::string> fPHist;
258  std::vector<std::string> fThetaXzYzHist;
259 
260  std::vector<std::unique_ptr<TH1>> hPHist;
261  std::vector<std::unique_ptr<TH2>>
263  // FYI - thetaxz and thetayz are related to standard polar angles as follows:
264  // thetaxz = atan2(math.sin(theta) * cos(phi), cos(theta))
265  // thetayz = asin(sin(theta) * sin(phi));
266 
267  CLHEP::HepRandomEngine& fEngine;
268 
270  static std::map<int, std::string> makeParticleSelectionModeNames();
271 
273  static std::map<int, std::string> makeDistributionNames();
274 
276  void beginRun(art::Run& run) override;
277 
279  void setup();
280 
303  template <typename OptionList>
304  static auto selectOption(std::string Option, OptionList const& allowedOptions)
305  -> decltype(auto);
306 
320  template <typename OptionList>
321  static std::string presentOptions(
322  OptionList const& allowedOptions,
323  bool printKey,
324  std::initializer_list<typename OptionList::value_type::first_type> exclude);
325 
326  template <typename OptionList>
327  static std::string presentOptions(OptionList const& allowedOptions, bool printKey = true)
328  {
329  return presentOptions(allowedOptions, printKey, {});
330  }
331 
333  template <typename OptionList>
334  static std::string optionName(typename OptionList::value_type::first_type optionKey,
335  OptionList const& allowedOptions,
336  std::string defName = "<unknown>");
337 
338  }; // class SingleGen
339 }
340 
341 namespace evgen {
342 
343  std::map<int, std::string> SingleGen::makeParticleSelectionModeNames()
344  {
345  std::map<int, std::string> names;
346  names[int(kSelectAllParts)] = "all";
347  names[int(kSelectOneRandPart)] = "singleRandom";
348  return names;
349  } // SingleGen::makeParticleSelectionModeNames()
350 
351  std::map<int, std::string> SingleGen::makeDistributionNames()
352  {
353  std::map<int, std::string> names;
354  names[int(kUNIF)] = "uniform";
355  names[int(kGAUS)] = "Gaussian";
356  names[int(kHIST)] = "histograms";
357  return names;
358  } // SingleGen::makeDistributionNames()
359 
360  const std::map<int, std::string> SingleGen::ParticleSelectionModeNames =
362  const std::map<int, std::string> SingleGen::DistributionNames =
364 
365  template <typename OptionList>
366  auto SingleGen::selectOption(std::string Option, OptionList const& allowedOptions)
367  -> decltype(auto)
368  {
369  using key_type = typename OptionList::value_type::first_type;
370  using tolower_type = int (*)(int);
371  auto toLower = [](auto const& S) {
372  std::string s;
373  s.reserve(S.size());
374  std::transform(S.cbegin(), S.cend(), std::back_inserter(s), (tolower_type)&std::tolower);
375  return s;
376  };
377  auto option = toLower(Option);
378  for (auto const& candidate : allowedOptions) {
379  if (toLower(candidate.second) == option) return candidate.first;
380  }
381  try {
382  std::size_t end;
383  key_type num = std::stoi(Option, &end);
384  if (allowedOptions.count(num) && (end == Option.length())) return num;
385  }
386  catch (std::invalid_argument const&) {
387  }
388  throw std::runtime_error("Option '" + Option + "' not supported.");
389  } // SingleGen::selectOption()
390 
391  template <typename OptionList>
393  OptionList const& allowedOptions,
394  bool printKey /* = true */,
395  std::initializer_list<typename OptionList::value_type::first_type> exclude /* = {} */
396  )
397  {
398  std::string msg;
399 
400  unsigned int n = 0;
401  for (auto const& option : allowedOptions) {
402  auto const& key = option.first;
403  if (std::find(exclude.begin(), exclude.end(), key) != exclude.end()) continue;
404  if (n++ > 0) msg += ", ";
405  msg += '\"' + std::string(option.second) + '\"';
406  if (printKey) msg += " (" + std::to_string(key) + ")";
407  } // for
408  return msg;
409  } // SingleGen::presentOptions()
410 
411  template <typename OptionList>
412  std::string SingleGen::optionName(typename OptionList::value_type::first_type optionKey,
413  OptionList const& allowedOptions,
414  std::string defName /* = "<unknown>" */
415  )
416  {
417  auto iOption = allowedOptions.find(optionKey);
418  return (iOption != allowedOptions.end()) ? iOption->second : defName;
419  } // SingleGen::optionName()
420 
421  //____________________________________________________________________________
422  bool SingleGen::Config::fromHistogram(std::string const& key) const
423  {
425  } // SingleGen::Config::fromHistogram()
426 
427  double SingleGen::particle_mass(int particle_id)
428  {
429  constexpr double kAlphaMass = 3.727379240;
430  if (particle_id == 1000020040) {
431  return kAlphaMass; // alpha particles
432  }
433 
434  static TDatabasePDG pdgt;
435  TParticlePDG* pdgp = pdgt.GetParticle(particle_id);
436  return pdgp ? pdgp->Mass() : 0.;
437  }
438 
439  simb::MCParticle SingleGen::mc_particle(int track_id, int particle_id, double mass)
440  {
441  if (particle_id == 1000020040) {
442  return simb::MCParticle{track_id, particle_id, "primary", -1, mass, 1};
443  }
444  return simb::MCParticle{track_id, particle_id, "primary"};
445  }
446 
447  //____________________________________________________________________________
449  : EDProducer{config}
450  , fMode(selectOption(config().ParticleSelectionMode(), ParticleSelectionModeNames))
451  , fPadOutVectors(config().PadOutVectors())
452  , fPDG(config().PDG())
453  , fP0(config().P0())
454  , fSigmaP(config().SigmaP())
455  , fPDist(selectOption(config().PDist(), DistributionNames))
456  , fX0(config().X0())
457  , fY0(config().Y0())
458  , fZ0(config().Z0())
459  , fT0(config().T0())
460  , fSigmaX(config().SigmaX())
461  , fSigmaY(config().SigmaY())
462  , fSigmaZ(config().SigmaZ())
463  , fSigmaT(config().SigmaT())
464  , fPosDist(selectOption(config().PosDist(), DistributionNames))
465  , fTDist(selectOption(config().TDist(), DistributionNames))
466  , fSingleVertex(config().SingleVertex())
467  , fTheta0XZ(config().Theta0XZ())
468  , fTheta0YZ(config().Theta0YZ())
469  , fSigmaThetaXZ(config().SigmaThetaXZ())
470  , fSigmaThetaYZ(config().SigmaThetaYZ())
471  , fAngleDist(selectOption(config().AngleDist(), DistributionNames))
472  , fHistFileName(config().HistogramFile())
473  , fPHist(config().PHist())
474  , fThetaXzYzHist(config().ThetaXzYzHist())
476  ->registerAndSeedEngine(createEngine(config().Seed()), "", "", config().Seed()))
477  {
478  setup();
479  mf::LogInfo("SingleGen") << "Seed set to " << fEngine.getSeed();
480 
481  produces<std::vector<simb::MCTruth>>();
482  produces<sumdata::RunData, art::InRun>();
483  }
484 
485  //____________________________________________________________________________
487  {
488  geo::GeometryCore const& geom = *(lar::providerFrom<geo::Geometry>());
489  run.put(std::make_unique<sumdata::RunData>(geom.DetectorName()), art::fullRun());
490  }
491 
492  //____________________________________________________________________________
494  {
495  // do not put seed in reconfigure because we don't want to reset
496  // the seed midstream
497  std::vector<std::string> vlist(15);
498  vlist[0] = "PDG";
499  vlist[1] = "P0";
500  vlist[2] = "SigmaP";
501  vlist[3] = "X0";
502  vlist[4] = "Y0";
503  vlist[5] = "Z0";
504  vlist[6] = "SigmaX";
505  vlist[7] = "SigmaY";
506  vlist[8] = "SigmaZ";
507  vlist[9] = "Theta0XZ";
508  vlist[10] = "Theta0YZ";
509  vlist[11] = "SigmaThetaXZ";
510  vlist[12] = "SigmaThetaYZ";
511  vlist[13] = "T0";
512  vlist[14] = "SigmaT";
513 
514  // vlist[15] = "PHist";
515  // vlist[16] = "ThetaHist";
516  // vlist[17] = "PhiHist";
517 
518  // begin tests for multiple particle error possibilities
519  std::string list;
520  if (fPDist != kHIST) {
521  if (!this->PadVector(fP0)) { list.append(vlist[1].append(", \n")); }
522  if (!this->PadVector(fSigmaP)) { list.append(vlist[2].append(", \n")); }
523  }
524  if (!this->PadVector(fX0)) { list.append(vlist[3].append(", \n")); }
525  if (!this->PadVector(fY0)) { list.append(vlist[4].append(", \n")); }
526  if (!this->PadVector(fZ0)) { list.append(vlist[5].append(", \n")); }
527  if (!this->PadVector(fSigmaX)) { list.append(vlist[6].append(", \n")); }
528  if (!this->PadVector(fSigmaY)) { list.append(vlist[7].append(", \n")); }
529  if (!this->PadVector(fSigmaZ)) { list.append(vlist[8].append(", \n")); }
530  if (!this->PadVector(fTheta0XZ)) { list.append(vlist[9].append(", \n")); }
531  if (!this->PadVector(fTheta0YZ)) { list.append(vlist[10].append(", \n")); }
532  if (!this->PadVector(fSigmaThetaXZ)) { list.append(vlist[11].append(", \n")); }
533  if (!this->PadVector(fSigmaThetaYZ)) { list.append(vlist[12].append(" \n")); }
534  if (!this->PadVector(fT0)) { list.append(vlist[13].append(", \n")); }
535  if (!this->PadVector(fSigmaT)) { list.append(vlist[14].append(", \n")); }
536 
537  if (list.size() > 0)
538  throw cet::exception("SingleGen")
539  << "The " << list << "\n vector(s) defined in the fhicl files has/have "
540  << "a different size than the PDG vector "
541  << "\n and it has (they have) more than one value, "
542  << "\n disallowing sensible padding "
543  << " and/or you have set fPadOutVectors to false. \n";
544 
545  if (fPDG.size() > 1 && fPadOutVectors) this->printVecs(vlist);
546 
547  // If needed, get histograms for momentum and angle distributions
548  TFile* histFile = nullptr;
549  if (!fHistFileName.empty()) {
550  if (fHistFileName[0] == '/') {
551  // We have an absolute path, use given name exactly.
552  if (cet::file_exists(fHistFileName)) {
553  histFile = new TFile(fHistFileName.c_str());
554  if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
555  delete histFile;
556  histFile = nullptr;
558  << "Cannot open ROOT file specified in parameter HistogramFile: \"" << fHistFileName
559  << "\"";
560  }
561  }
562  else {
564  << "ROOT file specified in parameter HistogramFile: \"" << fHistFileName
565  << "\" does not exist!";
566  }
567  }
568  else {
569  // We have a relative path, search starting from current directory.
570  std::string relative_filename{"./"};
571  relative_filename += fHistFileName;
572  if (cet::file_exists(relative_filename)) {
573  histFile = new TFile(relative_filename.c_str());
574  if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
575  delete histFile;
576  histFile = nullptr;
578  << "Cannot open ROOT file found using relative path and originally specified in "
579  "parameter HistogramFile: \""
580  << relative_filename << '"';
581  }
582  }
583  else {
584  cet::search_path sp{"FW_SEARCH_PATH"};
585  std::string found_filename;
586  auto found = sp.find_file(fHistFileName, found_filename);
587  if (!found) {
589  << "Cannot find ROOT file in current directory nor on FW_SEARCH_PATH specified in "
590  "parameter HistogramFile: \""
591  << fHistFileName << '"';
592  }
593  histFile = new TFile(found_filename.c_str());
594  if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
595  delete histFile;
596  histFile = nullptr;
598  << "Cannot open ROOT file found on FW_SEARCH_PATH and originally specified in "
599  "parameter HistogramFile: \""
600  << found_filename << '"';
601  }
602  }
603  }
604  }
605 
606  //
607  // deal with position distribution
608  //
609  switch (fPosDist) {
610  case kGAUS:
611  case kUNIF: break; // supported, no further action needed
612  default:
614  << "Position distribution of type '" << optionName(fPosDist, DistributionNames) << "' ("
615  << std::to_string(fPosDist) << ") is not supported.";
616  } // switch(fPosDist)
617 
618  //
619  // deal with time distribution
620  //
621  switch (fTDist) {
622  case kGAUS:
623  case kUNIF: break; // supported, no further action needed
624  default:
626  << "Time distribution of type '" << optionName(fTDist, DistributionNames) << "' ("
627  << std::to_string(fTDist) << ") is not supported.";
628  } // switch(fTDist)
629 
630  //
631  // deal with momentum distribution
632  //
633  switch (fPDist) {
634  case kHIST:
635  if (fPHist.size() != fPDG.size()) {
637  << fPHist.size() << " momentum histograms to describe " << fPDG.size()
638  << " particle types...";
639  }
640  hPHist.reserve(fPHist.size());
641  for (auto const& histName : fPHist) {
642  TH1* pHist = dynamic_cast<TH1*>(histFile->Get(histName.c_str()));
643  if (!pHist) {
645  << "Failed to read momentum histogram '" << histName << "' from '"
646  << histFile->GetPath() << "\'";
647  }
648  pHist->SetDirectory(nullptr); // make it independent of the input file
649  hPHist.emplace_back(pHist);
650  } // for
651  break;
652  default: // supported, no further action needed
653  break;
654  } // switch(fPDist)
655 
656  switch (fAngleDist) {
657  case kHIST:
658  if (fThetaXzYzHist.size() != fPDG.size()) {
660  << fThetaXzYzHist.size() << " direction histograms to describe " << fPDG.size()
661  << " particle types...";
662  }
663  hThetaXzYzHist.reserve(fThetaXzYzHist.size());
664  for (auto const& histName : fThetaXzYzHist) {
665  TH2* pHist = dynamic_cast<TH2*>(histFile->Get(histName.c_str()));
666  if (!pHist) {
668  << "Failed to read direction histogram '" << histName << "' from '"
669  << histFile->GetPath() << "\'";
670  }
671  pHist->SetDirectory(nullptr); // make it independent of the input file
672  hThetaXzYzHist.emplace_back(pHist);
673  } // for
674  default: // supported, no further action needed
675  break;
676  } // switch(fAngleDist)
677 
678  delete histFile;
679  }
680 
681  //____________________________________________________________________________
682  bool SingleGen::PadVector(std::vector<double>& vec)
683  {
684  // check if the vec has the same size as fPDG
685  if (vec.size() != fPDG.size()) {
686  // if not padding out the vectors always cause an
687  // exception to be thrown if the vector in question
688  // is not the same size as the fPDG vector
689  // the exception is thrown in the reconfigure method
690  // that calls this one
691  if (!fPadOutVectors)
692  return false;
693  else if (fPadOutVectors) {
694  // if padding of vectors is desired but the vector in
695  // question has more than one entry it isn't clear
696  // what the padded values should be so cause
697  // an exception
698  if (vec.size() != 1) return false;
699 
700  // pad it out
701  vec.resize(fPDG.size(), vec[0]);
702 
703  } // end if padding out vectors
704  } // end if the vector size is not the same as fPDG
705 
706  return true;
707  }
708 
709  //____________________________________________________________________________
711  {
712 
714  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
715 
716  simb::MCTruth truth;
718  Sample(truth);
719 
720  MF_LOG_DEBUG("SingleGen") << truth;
721 
722  truthcol->push_back(truth);
723 
724  evt.put(std::move(truthcol));
725 
726  return;
727  }
728 
729  //____________________________________________________________________________
730  // Draw the type, momentum and position of a single particle from the
731  // FCIHL description
732  void SingleGen::SampleOne(unsigned int i, simb::MCTruth& mct)
733  {
734 
735  CLHEP::RandFlat flat(fEngine);
736  CLHEP::RandGaussQ gauss(fEngine);
737 
738  // Choose momentum
739  double p = 0.0;
740  double m = 0.0;
741  if (fPDist == kGAUS) { p = gauss.fire(fP0[i], fSigmaP[i]); }
742  else if (fPDist == kHIST) {
743  p = SelectFromHist(*(hPHist[i]));
744  }
745  else { // if (fPDist == kUNIF) {
746  p = fP0[i] + fSigmaP[i] * (2.0 * flat.fire() - 1.0);
747  }
748  // else {std::cout << "do not understand the value of PDist!";}
749  m = particle_mass(fPDG[i]);
750 
751  // Choose position
752  TVector3 x;
753  if (fPosDist == kGAUS) {
754  x[0] = gauss.fire(fX0[i], fSigmaX[i]);
755  ;
756  x[1] = gauss.fire(fY0[i], fSigmaY[i]);
757  x[2] = gauss.fire(fZ0[i], fSigmaZ[i]);
758  }
759  else {
760  x[0] = fX0[i] + fSigmaX[i] * (2.0 * flat.fire() - 1.0);
761  x[1] = fY0[i] + fSigmaY[i] * (2.0 * flat.fire() - 1.0);
762  x[2] = fZ0[i] + fSigmaZ[i] * (2.0 * flat.fire() - 1.0);
763  }
764 
765  double t = 0.;
766  if (fTDist == kGAUS) { t = gauss.fire(fT0[i], fSigmaT[i]); }
767  else {
768  t = fT0[i] + fSigmaT[i] * (2.0 * flat.fire() - 1.0);
769  }
770 
771  TLorentzVector pos(x[0], x[1], x[2], t);
772 
773  // Choose angles
774  double thxz = 0;
775  double thyz = 0;
776 
777  double thyzrads = 0;
778  double thyzradsplussigma = 0;
779  double thyzradsminussigma = 0;
780 
781  if (fAngleDist == kGAUS) {
782  thxz = gauss.fire(fTheta0XZ[i], fSigmaThetaXZ[i]);
783  thyz = gauss.fire(fTheta0YZ[i], fSigmaThetaYZ[i]);
784  }
785  else if (fAngleDist == kHIST) { // Select thetaxz and thetayz from histogram
786  double thetaxz = 0;
787  double thetayz = 0;
788  SelectFromHist(*(hThetaXzYzHist[i]), thetaxz, thetayz);
789  thxz = (180. / M_PI) * thetaxz;
790  thyz = (180. / M_PI) * thetayz;
791  }
792  else {
793 
794  // Choose angles flat in phase space, which is flat in theta_xz
795  // and flat in sin(theta_yz).
796 
797  thxz = fTheta0XZ[i] + fSigmaThetaXZ[i] * (2.0 * flat.fire() - 1.0);
798 
799  thyzrads = std::asin(std::sin(
800  (M_PI / 180.) *
801  (fTheta0YZ
802  [i]))); //Taking asin of sin gives value between -Pi/2 and Pi/2 regardless of user input
803  thyzradsplussigma =
804  TMath::Min((thyzrads + ((M_PI / 180.) * fabs(fSigmaThetaYZ[i]))), M_PI / 2.);
805  thyzradsminussigma =
806  TMath::Max((thyzrads - ((M_PI / 180.) * fabs(fSigmaThetaYZ[i]))), -M_PI / 2.);
807 
808  //uncomment line to print angular variation info
809  //std::cout << "Central angle: " << (180./M_PI)*thyzrads << " Max angle: " << (180./M_PI)*thyzradsplussigma << " Min angle: " << (180./M_PI)*thyzradsminussigma << std::endl;
810 
811  double sinthyzmin = std::sin(thyzradsminussigma);
812  double sinthyzmax = std::sin(thyzradsplussigma);
813  double sinthyz = sinthyzmin + flat.fire() * (sinthyzmax - sinthyzmin);
814  thyz = (180. / M_PI) * std::asin(sinthyz);
815  }
816 
817  double thxzrad = thxz * M_PI / 180.0;
818  double thyzrad = thyz * M_PI / 180.0;
819 
820  TLorentzVector pvec(p * std::cos(thyzrad) * std::sin(thxzrad),
821  p * std::sin(thyzrad),
822  p * std::cos(thxzrad) * std::cos(thyzrad),
823  std::sqrt(p * p + m * m));
824 
825  // set track id to -i as these are all primary particles and have id <= 0
826  int trackid = -1 * (i + 1);
827  std::string primary("primary");
828 
829  auto part = mc_particle(trackid, fPDG[i], m);
830  part.AddTrajectoryPoint(pos, pvec);
831  mct.Add(part);
832 
833  //std::cout << "Px: " << pvec.Px() << " Py: " << pvec.Py() << " Pz: " << pvec.Pz() << std::endl;
834  //std::cout << "x: " << pos.X() << " y: " << pos.Y() << " z: " << pos.Z() << " time: " << pos.T() << std::endl;
835  //std::cout << "YZ Angle: " << (thyzrad * (180./M_PI)) << " XZ Angle: " << (thxzrad * (180./M_PI)) << std::endl;
836  }
837 
838  //____________________________________________________________________________
839  // Draw the type, momentum and position for all particles from the
840  // FCIHL description. Start positions will all match but momenta and angles drawn from
841  // distributions defined in the fhicls
843  {
844 
845  CLHEP::RandFlat flat(fEngine);
846  CLHEP::RandGaussQ gauss(fEngine);
847 
848  // Choose position
849  TVector3 x;
850  if (fPosDist == kGAUS) {
851  x[0] = gauss.fire(fX0[0], fSigmaX[0]);
852  ;
853  x[1] = gauss.fire(fY0[0], fSigmaY[0]);
854  x[2] = gauss.fire(fZ0[0], fSigmaZ[0]);
855  }
856  else {
857  x[0] = fX0[0] + fSigmaX[0] * (2.0 * flat.fire() - 1.0);
858  x[1] = fY0[0] + fSigmaY[0] * (2.0 * flat.fire() - 1.0);
859  x[2] = fZ0[0] + fSigmaZ[0] * (2.0 * flat.fire() - 1.0);
860  }
861 
862  double t = 0.;
863  if (fTDist == kGAUS) { t = gauss.fire(fT0[0], fSigmaT[0]); }
864  else {
865  t = fT0[0] + fSigmaT[0] * (2.0 * flat.fire() - 1.0);
866  }
867 
868  TLorentzVector pos(x[0], x[1], x[2], t);
869 
870  // loop through particles and select momenta and angles
871  for (unsigned int i(0); i < fPDG.size(); ++i) {
872  // Choose momentum
873  double p = 0.0;
874  double m = 0.0;
875  if (fPDist == kGAUS) { p = gauss.fire(fP0[i], fSigmaP[i]); }
876  else if (fPDist == kHIST) {
877  p = SelectFromHist(*(hPHist[i]));
878  }
879  else {
880  p = fP0[i] + fSigmaP[i] * (2.0 * flat.fire() - 1.0);
881  }
882 
883  m = particle_mass(fPDG[i]);
884 
885  // Choose angles
886  double thxz = 0;
887  double thyz = 0;
888 
889  double thyzrads = 0;
890  double thyzradsplussigma = 0;
891  double thyzradsminussigma = 0;
892 
893  if (fAngleDist == kGAUS) {
894  thxz = gauss.fire(fTheta0XZ[i], fSigmaThetaXZ[i]);
895  thyz = gauss.fire(fTheta0YZ[i], fSigmaThetaYZ[i]);
896  }
897  else if (fAngleDist == kHIST) {
898  double thetaxz = 0;
899  double thetayz = 0;
900  SelectFromHist(*(hThetaXzYzHist[i]), thetaxz, thetayz);
901  thxz = (180. / M_PI) * thetaxz;
902  thyz = (180. / M_PI) * thetayz;
903  }
904  else {
905 
906  // Choose angles flat in phase space, which is flat in theta_xz
907  // and flat in sin(theta_yz).
908 
909  thxz = fTheta0XZ[i] + fSigmaThetaXZ[i] * (2.0 * flat.fire() - 1.0);
910 
911  thyzrads = std::asin(std::sin(
912  (M_PI / 180.) *
913  (fTheta0YZ
914  [i]))); //Taking asin of sin gives value between -Pi/2 and Pi/2 regardless of user input
915  thyzradsplussigma =
916  TMath::Min((thyzrads + ((M_PI / 180.) * fabs(fSigmaThetaYZ[i]))), M_PI / 2.);
917  thyzradsminussigma =
918  TMath::Max((thyzrads - ((M_PI / 180.) * fabs(fSigmaThetaYZ[i]))), -M_PI / 2.);
919 
920  //uncomment line to print angular variation info
921  //std::cout << "Central angle: " << (180./M_PI)*thyzrads << " Max angle: " << (180./M_PI)*thyzradsplussigma << " Min angle: " << (180./M_PI)*thyzradsminussigma << std::endl;
922 
923  double sinthyzmin = std::sin(thyzradsminussigma);
924  double sinthyzmax = std::sin(thyzradsplussigma);
925  double sinthyz = sinthyzmin + flat.fire() * (sinthyzmax - sinthyzmin);
926  thyz = (180. / M_PI) * std::asin(sinthyz);
927  }
928 
929  double thxzrad = thxz * M_PI / 180.0;
930  double thyzrad = thyz * M_PI / 180.0;
931 
932  TLorentzVector pvec(p * std::cos(thyzrad) * std::sin(thxzrad),
933  p * std::sin(thyzrad),
934  p * std::cos(thxzrad) * std::cos(thyzrad),
935  std::sqrt(p * p + m * m));
936 
937  // set track id to -i as these are all primary particles and have id <= 0
938  int trackid = -1 * (i + 1);
939  std::string primary("primary");
940 
941  auto part = mc_particle(trackid, fPDG[i], m);
942  part.AddTrajectoryPoint(pos, pvec);
943  mct.Add(part);
944 
945  //std::cout << "Px: " << pvec.Px() << " Py: " << pvec.Py() << " Pz: " << pvec.Pz() << std::endl;
946  //std::cout << "x: " << pos.X() << " y: " << pos.Y() << " z: " << pos.Z() << " time: " << pos.T() << std::endl;
947  //std::cout << "YZ Angle: " << (thyzrad * (180./M_PI)) << " XZ Angle: " << (thxzrad * (180./M_PI)) << std::endl;
948  }
949  }
950 
951  //____________________________________________________________________________
953  {
954 
955  switch (fMode) {
956  case 0: // List generation mode: every event will have one of each
957  // particle species in the fPDG array
958  if (fSingleVertex) { SampleMany(mct); }
959  else {
960  for (unsigned int i = 0; i < fPDG.size(); ++i) {
961  SampleOne(i, mct);
962  } //end loop over particles
963  }
964  break;
965  case 1: // Random selection mode: every event will exactly one particle
966  // selected randomly from the fPDG array
967  {
968  CLHEP::RandFlat flat(fEngine);
969 
970  unsigned int i = flat.fireInt(fPDG.size());
971  SampleOne(i, mct);
972  } break;
973  default:
974  mf::LogWarning("UnrecognizeOption")
975  << "SingleGen does not recognize ParticleSelectionMode " << fMode;
976  break;
977  } // switch on fMode
978 
979  return;
980  }
981 
982  //____________________________________________________________________________
983  void SingleGen::printVecs(std::vector<std::string> const& list)
984  {
985 
986  mf::LogInfo("SingleGen") << " You are using vector values for SingleGen configuration.\n "
987  << " Some of the configuration vectors may have been padded out ,"
988  << " because they (weren't) as long as the pdg vector"
989  << " in your configuration. \n"
990  << " The new input particle configuration is:\n";
991 
992  std::string values;
993  for (size_t i = 0; i <= 1; ++i) { // list.size(); ++i){
994 
995  values.append(list[i]);
996  values.append(": [ ");
997 
998  for (size_t e = 0; e < fPDG.size(); ++e) {
999  std::stringstream buf;
1000  buf.width(10);
1001  if (i == 0) buf << fPDG[e] << ", ";
1002  buf.precision(5);
1003  if (i == 1) buf << fP0[e] << ", ";
1004  if (i == 2) buf << fSigmaP[e] << ", ";
1005  if (i == 3) buf << fX0[e] << ", ";
1006  if (i == 4) buf << fY0[e] << ", ";
1007  if (i == 5) buf << fZ0[e] << ", ";
1008  if (i == 6) buf << fSigmaX[e] << ", ";
1009  if (i == 7) buf << fSigmaY[e] << ", ";
1010  if (i == 8) buf << fSigmaZ[e] << ", ";
1011  if (i == 9) buf << fTheta0XZ[e] << ", ";
1012  if (i == 10) buf << fTheta0YZ[e] << ", ";
1013  if (i == 11) buf << fSigmaThetaXZ[e] << ", ";
1014  if (i == 12) buf << fSigmaThetaYZ[e] << ", ";
1015  if (i == 13) buf << fT0[e] << ", ";
1016  if (i == 14) buf << fSigmaT[e] << ", ";
1017  values.append(buf.str());
1018  }
1019 
1020  values.erase(values.find_last_of(","));
1021  values.append(" ] \n");
1022 
1023  } // end loop over vector names in list
1024 
1025  mf::LogInfo("SingleGen") << values;
1026 
1027  return;
1028  }
1029 
1030  //____________________________________________________________________________
1031  double SingleGen::SelectFromHist(const TH1& h) // select from a 1D histogram
1032  {
1033  CLHEP::RandFlat flat(fEngine);
1034 
1035  double throw_value = h.Integral() * flat.fire();
1036  double cum_value(0);
1037  for (int i(0); i < h.GetNbinsX() + 1; ++i) {
1038  cum_value += h.GetBinContent(i);
1039  if (throw_value < cum_value) { return flat.fire() * h.GetBinWidth(i) + h.GetBinLowEdge(i); }
1040  }
1041  return throw_value; // for some reason we've gone through all bins and failed?
1042  }
1043  //____________________________________________________________________________
1044  void SingleGen::SelectFromHist(const TH2& h, double& x, double& y) // select from a 2D histogram
1045  {
1046  CLHEP::RandFlat flat(fEngine);
1047 
1048  double throw_value = h.Integral() * flat.fire();
1049  double cum_value(0);
1050  for (int i(0); i < h.GetNbinsX() + 1; ++i) {
1051  for (int j(0); j < h.GetNbinsY() + 1; ++j) {
1052  cum_value += h.GetBinContent(i, j);
1053  if (throw_value < cum_value) {
1054  x = flat.fire() * h.GetXaxis()->GetBinWidth(i) + h.GetXaxis()->GetBinLowEdge(i);
1055  y = flat.fire() * h.GetYaxis()->GetBinWidth(j) + h.GetYaxis()->GetBinLowEdge(j);
1056  return;
1057  }
1058  }
1059  }
1060  return; // for some reason we've gone through all bins and failed?
1061  }
1062  //____________________________________________________________________________
1063 
1064 } //end namespace evgen
1065 
1066 namespace evgen {
1067 
1069 
1070 } //end namespace evgen
Float_t x
Definition: compare.C:6
int fPDist
How to distribute momenta (gaus or uniform)
static const std::map< int, std::string > DistributionNames
Names of all distribution modes.
base_engine_t & createEngine(seed_t seed)
int fTDist
How to distribute t (gaus, or uniform)
module to produce single or multiple specified particles in the detector
Utilities related to art service access.
bool PadVector(std::vector< double > &vec)
static std::string presentOptions(OptionList const &allowedOptions, bool printKey, std::initializer_list< typename OptionList::value_type::first_type > exclude)
Returns a string describing all options in the list.
constexpr auto fullRun()
fhicl::Sequence< double > SigmaT
static std::map< int, std::string > makeParticleSelectionModeNames()
Returns a vector with the name of particle selection mode keywords.
static std::map< int, std::string > makeDistributionNames()
Returns a vector with the name of distribution keywords.
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
fhicl::Sequence< double > SigmaP
std::vector< std::unique_ptr< TH2 > > hThetaXzYzHist
actual TH1 for momentum distributions
std::vector< double > fSigmaT
Variation in t position (ns)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
fhicl::Sequence< double > SigmaThetaXZ
Float_t y
Definition: compare.C:6
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
std::vector< std::unique_ptr< TH1 > > hPHist
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Run.h:121
static std::string presentOptions(OptionList const &allowedOptions, bool printKey=true)
std::vector< std::string > fThetaXzYzHist
name of histogram for thetaxz/thetayz distribution
fhicl::Sequence< double > Y0
void SampleOne(unsigned int i, simb::MCTruth &mct)
bool fSingleVertex
if true - all particles produced at the same location
fhicl::Sequence< double > SigmaThetaYZ
fhicl::Sequence< double > SigmaZ
Particle class.
fhicl::Sequence< std::string > ThetaXzYzHist
int fPosDist
How to distribute xyz (gaus, or uniform)
void printVecs(std::vector< std::string > const &list)
fhicl::Sequence< double > P0
Definition: Run.h:37
fhicl::Sequence< double > SigmaY
std::vector< double > fSigmaZ
Variation in z position (cm)
fhicl::Sequence< double > X0
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
static constexpr int kSelectOneRandPart
fhicl::Sequence< double > Theta0YZ
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
Access the description of the physical detector geometry.
fhicl::Sequence< double > T0
fhicl::Atom< std::string > HistogramFile
std::vector< double > fT0
Central t position (ns) in world coordinates.
TString part[npart]
Definition: Style.C:32
fhicl::Atom< std::string > AngleDist
void produce(art::Event &evt) override
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
decltype(auto) constexpr to_string(T &&obj)
ADL-aware version of std::to_string.
std::vector< double > fP0
Central momentum (GeV/c) to generate.
std::vector< double > fTheta0YZ
Angle in YZ plane (degrees)
single particles thrown at the detector
Definition: MCTruth.h:26
decltype(auto) values(Coll &&coll)
Range-for loop helper iterating across the values of the specified collection.
fhicl::Atom< bool > SingleVertex
fhicl::Sequence< double > SigmaX
static constexpr seed_t InvalidSeed
An invalid seed.
fhicl::Atom< std::string > ParticleSelectionMode
static constexpr int kSelectAllParts
One particle per entry is generated.
double SelectFromHist(const TH1 &h)
void Sample(simb::MCTruth &mct)
fhicl::Sequence< int > PDG
Description of the physical geometry of one entire detector.
Definition: GeometryCore.h:91
static auto selectOption(std::string Option, OptionList const &allowedOptions) -> decltype(auto)
Parses an option string and returns the corresponding option number.
fhicl::Sequence< double > Theta0XZ
fhicl::Sequence< double > Z0
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
std::vector< double > fY0
Central y position (cm) in world coordinates.
std::vector< int > fPDG
PDG code of particles to generate.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::vector< double > fSigmaThetaYZ
Variation in angle in YZ plane.
std::vector< double > fSigmaP
Variation in momenta (GeV/c)
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
fhicl::Sequence< std::string > PHist
std::vector< double > fZ0
Central z position (cm) in world coordinates.
std::vector< std::string > fPHist
name of histogram of momenta
std::vector< double > fSigmaX
Variation in x position (cm)
fhicl::Atom< std::string > PosDist
static constexpr int kGAUS
Gaussian distribution.
bool fromHistogram(std::string const &key) const
Returns whether the specified mode is an histogram distribution.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double particle_mass(int particle_id)
Returns the particle mass for the specified PDG code.
#define MF_LOG_DEBUG(id)
std::vector< double > fX0
Central x position (cm) in world coordinates.
static std::string optionName(typename OptionList::value_type::first_type optionKey, OptionList const &allowedOptions, std::string defName="<unknown>")
Returns the name of the specified option key, or defName if not known.
Char_t n[5]
static const std::map< int, std::string > ParticleSelectionModeNames
Names of all particle selection modes.
std::string fHistFileName
Filename containing histogram of momenta.
int fAngleDist
How to distribute angles (gaus, uniform)
SingleGen(Parameters const &config)
static constexpr int kUNIF
Uniform distribution.
std::vector< double > fTheta0XZ
Angle in XZ plane (degrees)
std::string const & DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:140
TCEvent evt
Definition: DataStructs.cxx:8
Event generator information.
Definition: MCTruth.h:32
void beginRun(art::Run &run) override
Act on begin of run: write "RunData" information (sumdata::RunData).
Float_t e
Definition: plot.C:35
static constexpr int kHIST
void setup()
Performs checks and initialization based on the current configuration.
fhicl::Atom< std::string > PDist
Event Generation using GENIE, cosmics or single particles.
fhicl::Atom< bool > PadOutVectors
simb::MCParticle mc_particle(int track_id, int particle_id, double mass)
Returns a simb::MCParticle object for the specified track and particle IDs.
fhicl::Atom< rndm::NuRandomService::seed_t > Seed
std::vector< double > fSigmaThetaXZ
Variation in angle in XZ plane.
void SampleMany(simb::MCTruth &mct)
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
fhicl::Atom< std::string > TDist
CLHEP::HepRandomEngine & fEngine
actual TH2 for angle distributions - Xz on x axis .
std::vector< double > fSigmaY
Variation in y position (cm)