LArSoft  v09_90_00
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 (s) [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 (s) [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")};
180 
181  private:
183  bool fromHistogram(std::string const& key) const;
184 
185  }; // struct Config
186 
188 
189  explicit SingleGen(Parameters const& config);
190 
191  // This is called for each event.
192  void produce(art::Event& evt) override;
193 
194  private:
196  static const std::map<int, std::string> ParticleSelectionModeNames;
198  static const std::map<int, std::string> DistributionNames;
199 
200  void SampleOne(unsigned int i, simb::MCTruth& mct);
201  void SampleMany(simb::MCTruth& mct);
202  void Sample(simb::MCTruth& mct);
203  void printVecs(std::vector<std::string> const& list);
204  bool PadVector(std::vector<double>& vec);
205  double SelectFromHist(const TH1& h);
206  void SelectFromHist(const TH2& h, double& x, double& y);
207 
210 
211  static constexpr int kSelectAllParts = 0;
212  static constexpr int kSelectOneRandPart =
213  1;
214 
218 
219  static constexpr int kUNIF = 0;
220  static constexpr int kGAUS = 1;
221  static constexpr int kHIST = 2;
222 
224  int fMode;
225  bool fPadOutVectors;
228  std::vector<int> fPDG;
232  std::vector<double> fP0;
233  std::vector<double> fSigmaP;
234  int fPDist;
235  std::vector<double> fX0;
236  std::vector<double> fY0;
237  std::vector<double> fZ0;
238  std::vector<double> fT0;
239  std::vector<double> fSigmaX;
240  std::vector<double> fSigmaY;
241  std::vector<double> fSigmaZ;
242  std::vector<double> fSigmaT;
243  int fPosDist;
244  int fTDist;
246  std::vector<double> fTheta0XZ;
247  std::vector<double> fTheta0YZ;
248  std::vector<double> fSigmaThetaXZ;
249  std::vector<double> fSigmaThetaYZ;
251  std::string fHistFileName;
252  std::vector<std::string> fPHist;
253  std::vector<std::string> fThetaXzYzHist;
254 
255  std::vector<std::unique_ptr<TH1>> hPHist;
256  std::vector<std::unique_ptr<TH2>>
258  // FYI - thetaxz and thetayz are related to standard polar angles as follows:
259  // thetaxz = atan2(math.sin(theta) * cos(phi), cos(theta))
260  // thetayz = asin(sin(theta) * sin(phi));
261 
262  CLHEP::HepRandomEngine& fEngine;
263 
265  static std::map<int, std::string> makeParticleSelectionModeNames();
266 
268  static std::map<int, std::string> makeDistributionNames();
269 
271  void beginRun(art::Run& run) override;
272 
274  void setup();
275 
298  template <typename OptionList>
299  static auto selectOption(std::string Option, OptionList const& allowedOptions)
300  -> decltype(auto);
301 
315  template <typename OptionList>
316  static std::string presentOptions(
317  OptionList const& allowedOptions,
318  bool printKey,
319  std::initializer_list<typename OptionList::value_type::first_type> exclude);
320 
321  template <typename OptionList>
322  static std::string presentOptions(OptionList const& allowedOptions, bool printKey = true)
323  {
324  return presentOptions(allowedOptions, printKey, {});
325  }
326 
328  template <typename OptionList>
329  static std::string optionName(typename OptionList::value_type::first_type optionKey,
330  OptionList const& allowedOptions,
331  std::string defName = "<unknown>");
332 
333  }; // class SingleGen
334 }
335 
336 namespace evgen {
337 
338  std::map<int, std::string> SingleGen::makeParticleSelectionModeNames()
339  {
340  std::map<int, std::string> names;
341  names[int(kSelectAllParts)] = "all";
342  names[int(kSelectOneRandPart)] = "singleRandom";
343  return names;
344  } // SingleGen::makeParticleSelectionModeNames()
345 
346  std::map<int, std::string> SingleGen::makeDistributionNames()
347  {
348  std::map<int, std::string> names;
349  names[int(kUNIF)] = "uniform";
350  names[int(kGAUS)] = "Gaussian";
351  names[int(kHIST)] = "histograms";
352  return names;
353  } // SingleGen::makeDistributionNames()
354 
355  const std::map<int, std::string> SingleGen::ParticleSelectionModeNames =
357  const std::map<int, std::string> SingleGen::DistributionNames =
359 
360  template <typename OptionList>
361  auto SingleGen::selectOption(std::string Option, OptionList const& allowedOptions)
362  -> decltype(auto)
363  {
364  using key_type = typename OptionList::value_type::first_type;
365  using tolower_type = int (*)(int);
366  auto toLower = [](auto const& S) {
367  std::string s;
368  s.reserve(S.size());
369  std::transform(S.cbegin(), S.cend(), std::back_inserter(s), (tolower_type)&std::tolower);
370  return s;
371  };
372  auto option = toLower(Option);
373  for (auto const& candidate : allowedOptions) {
374  if (toLower(candidate.second) == option) return candidate.first;
375  }
376  try {
377  std::size_t end;
378  key_type num = std::stoi(Option, &end);
379  if (allowedOptions.count(num) && (end == Option.length())) return num;
380  }
381  catch (std::invalid_argument const&) {
382  }
383  throw std::runtime_error("Option '" + Option + "' not supported.");
384  } // SingleGen::selectOption()
385 
386  template <typename OptionList>
388  OptionList const& allowedOptions,
389  bool printKey /* = true */,
390  std::initializer_list<typename OptionList::value_type::first_type> exclude /* = {} */
391  )
392  {
393  std::string msg;
394 
395  unsigned int n = 0;
396  for (auto const& option : allowedOptions) {
397  auto const& key = option.first;
398  if (std::find(exclude.begin(), exclude.end(), key) != exclude.end()) continue;
399  if (n++ > 0) msg += ", ";
400  msg += '\"' + std::string(option.second) + '\"';
401  if (printKey) msg += " (" + std::to_string(key) + ")";
402  } // for
403  return msg;
404  } // SingleGen::presentOptions()
405 
406  template <typename OptionList>
407  std::string SingleGen::optionName(typename OptionList::value_type::first_type optionKey,
408  OptionList const& allowedOptions,
409  std::string defName /* = "<unknown>" */
410  )
411  {
412  auto iOption = allowedOptions.find(optionKey);
413  return (iOption != allowedOptions.end()) ? iOption->second : defName;
414  } // SingleGen::optionName()
415 
416  //____________________________________________________________________________
417  bool SingleGen::Config::fromHistogram(std::string const& key) const
418  {
420  } // SingleGen::Config::fromHistogram()
421 
422  //____________________________________________________________________________
424  : EDProducer{config}
425  , fMode(selectOption(config().ParticleSelectionMode(), ParticleSelectionModeNames))
426  , fPadOutVectors(config().PadOutVectors())
427  , fPDG(config().PDG())
428  , fP0(config().P0())
429  , fSigmaP(config().SigmaP())
430  , fPDist(selectOption(config().PDist(), DistributionNames))
431  , fX0(config().X0())
432  , fY0(config().Y0())
433  , fZ0(config().Z0())
434  , fT0(config().T0())
435  , fSigmaX(config().SigmaX())
436  , fSigmaY(config().SigmaY())
437  , fSigmaZ(config().SigmaZ())
438  , fSigmaT(config().SigmaT())
439  , fPosDist(selectOption(config().PosDist(), DistributionNames))
440  , fTDist(selectOption(config().TDist(), DistributionNames))
441  , fSingleVertex(config().SingleVertex())
442  , fTheta0XZ(config().Theta0XZ())
443  , fTheta0YZ(config().Theta0YZ())
444  , fSigmaThetaXZ(config().SigmaThetaXZ())
445  , fSigmaThetaYZ(config().SigmaThetaYZ())
446  , fAngleDist(selectOption(config().AngleDist(), DistributionNames))
447  , fHistFileName(config().HistogramFile())
448  , fPHist(config().PHist())
449  , fThetaXzYzHist(config().ThetaXzYzHist())
450  , fEngine(createEngine(0))
451  {
452  setup();
454  if (config().Seed(seed)) { fEngine.setSeed(seed, 0 /* dummy? */); }
455 
456  produces<std::vector<simb::MCTruth>>();
457  produces<sumdata::RunData, art::InRun>();
458  }
459 
460  //____________________________________________________________________________
462  {
463  geo::GeometryCore const& geom = *(lar::providerFrom<geo::Geometry>());
464  run.put(std::make_unique<sumdata::RunData>(geom.DetectorName()), art::fullRun());
465  }
466 
467  //____________________________________________________________________________
469  {
470  // do not put seed in reconfigure because we don't want to reset
471  // the seed midstream
472  std::vector<std::string> vlist(15);
473  vlist[0] = "PDG";
474  vlist[1] = "P0";
475  vlist[2] = "SigmaP";
476  vlist[3] = "X0";
477  vlist[4] = "Y0";
478  vlist[5] = "Z0";
479  vlist[6] = "SigmaX";
480  vlist[7] = "SigmaY";
481  vlist[8] = "SigmaZ";
482  vlist[9] = "Theta0XZ";
483  vlist[10] = "Theta0YZ";
484  vlist[11] = "SigmaThetaXZ";
485  vlist[12] = "SigmaThetaYZ";
486  vlist[13] = "T0";
487  vlist[14] = "SigmaT";
488 
489  // vlist[15] = "PHist";
490  // vlist[16] = "ThetaHist";
491  // vlist[17] = "PhiHist";
492 
493  // begin tests for multiple particle error possibilities
494  std::string list;
495  if (fPDist != kHIST) {
496  if (!this->PadVector(fP0)) { list.append(vlist[1].append(", \n")); }
497  if (!this->PadVector(fSigmaP)) { list.append(vlist[2].append(", \n")); }
498  }
499  if (!this->PadVector(fX0)) { list.append(vlist[3].append(", \n")); }
500  if (!this->PadVector(fY0)) { list.append(vlist[4].append(", \n")); }
501  if (!this->PadVector(fZ0)) { list.append(vlist[5].append(", \n")); }
502  if (!this->PadVector(fSigmaX)) { list.append(vlist[6].append(", \n")); }
503  if (!this->PadVector(fSigmaY)) { list.append(vlist[7].append(", \n")); }
504  if (!this->PadVector(fSigmaZ)) { list.append(vlist[8].append(", \n")); }
505  if (!this->PadVector(fTheta0XZ)) { list.append(vlist[9].append(", \n")); }
506  if (!this->PadVector(fTheta0YZ)) { list.append(vlist[10].append(", \n")); }
507  if (!this->PadVector(fSigmaThetaXZ)) { list.append(vlist[11].append(", \n")); }
508  if (!this->PadVector(fSigmaThetaYZ)) { list.append(vlist[12].append(" \n")); }
509  if (!this->PadVector(fT0)) { list.append(vlist[13].append(", \n")); }
510  if (!this->PadVector(fSigmaT)) { list.append(vlist[14].append(", \n")); }
511 
512  if (list.size() > 0)
513  throw cet::exception("SingleGen")
514  << "The " << list << "\n vector(s) defined in the fhicl files has/have "
515  << "a different size than the PDG vector "
516  << "\n and it has (they have) more than one value, "
517  << "\n disallowing sensible padding "
518  << " and/or you have set fPadOutVectors to false. \n";
519 
520  if (fPDG.size() > 1 && fPadOutVectors) this->printVecs(vlist);
521 
522  // If needed, get histograms for momentum and angle distributions
523  TFile* histFile = nullptr;
524  if (!fHistFileName.empty()) {
525  if (fHistFileName[0] == '/') {
526  // We have an absolute path, use given name exactly.
527  if (cet::file_exists(fHistFileName)) {
528  histFile = new TFile(fHistFileName.c_str());
529  if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
530  delete histFile;
531  histFile = nullptr;
533  << "Cannot open ROOT file specified in parameter HistogramFile: \"" << fHistFileName
534  << "\"";
535  }
536  }
537  else {
539  << "ROOT file specified in parameter HistogramFile: \"" << fHistFileName
540  << "\" does not exist!";
541  }
542  }
543  else {
544  // We have a relative path, search starting from current directory.
545  std::string relative_filename{"./"};
546  relative_filename += fHistFileName;
547  if (cet::file_exists(relative_filename)) {
548  histFile = new TFile(relative_filename.c_str());
549  if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
550  delete histFile;
551  histFile = nullptr;
553  << "Cannot open ROOT file found using relative path and originally specified in "
554  "parameter HistogramFile: \""
555  << relative_filename << '"';
556  }
557  }
558  else {
559  cet::search_path sp{"FW_SEARCH_PATH"};
560  std::string found_filename;
561  auto found = sp.find_file(fHistFileName, found_filename);
562  if (!found) {
564  << "Cannot find ROOT file in current directory nor on FW_SEARCH_PATH specified in "
565  "parameter HistogramFile: \""
566  << fHistFileName << '"';
567  }
568  histFile = new TFile(found_filename.c_str());
569  if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
570  delete histFile;
571  histFile = nullptr;
573  << "Cannot open ROOT file found on FW_SEARCH_PATH and originally specified in "
574  "parameter HistogramFile: \""
575  << found_filename << '"';
576  }
577  }
578  }
579  }
580 
581  //
582  // deal with position distribution
583  //
584  switch (fPosDist) {
585  case kGAUS:
586  case kUNIF: break; // supported, no further action needed
587  default:
589  << "Position distribution of type '" << optionName(fPosDist, DistributionNames) << "' ("
590  << std::to_string(fPosDist) << ") is not supported.";
591  } // switch(fPosDist)
592 
593  //
594  // deal with time distribution
595  //
596  switch (fTDist) {
597  case kGAUS:
598  case kUNIF: break; // supported, no further action needed
599  default:
601  << "Time distribution of type '" << optionName(fTDist, DistributionNames) << "' ("
602  << std::to_string(fTDist) << ") is not supported.";
603  } // switch(fTDist)
604 
605  //
606  // deal with momentum distribution
607  //
608  switch (fPDist) {
609  case kHIST:
610  if (fPHist.size() != fPDG.size()) {
612  << fPHist.size() << " momentum histograms to describe " << fPDG.size()
613  << " particle types...";
614  }
615  hPHist.reserve(fPHist.size());
616  for (auto const& histName : fPHist) {
617  TH1* pHist = dynamic_cast<TH1*>(histFile->Get(histName.c_str()));
618  if (!pHist) {
620  << "Failed to read momentum histogram '" << histName << "' from '"
621  << histFile->GetPath() << "\'";
622  }
623  pHist->SetDirectory(nullptr); // make it independent of the input file
624  hPHist.emplace_back(pHist);
625  } // for
626  break;
627  default: // supported, no further action needed
628  break;
629  } // switch(fPDist)
630 
631  switch (fAngleDist) {
632  case kHIST:
633  if (fThetaXzYzHist.size() != fPDG.size()) {
635  << fThetaXzYzHist.size() << " direction histograms to describe " << fPDG.size()
636  << " particle types...";
637  }
638  hThetaXzYzHist.reserve(fThetaXzYzHist.size());
639  for (auto const& histName : fThetaXzYzHist) {
640  TH2* pHist = dynamic_cast<TH2*>(histFile->Get(histName.c_str()));
641  if (!pHist) {
643  << "Failed to read direction histogram '" << histName << "' from '"
644  << histFile->GetPath() << "\'";
645  }
646  pHist->SetDirectory(nullptr); // make it independent of the input file
647  hThetaXzYzHist.emplace_back(pHist);
648  } // for
649  default: // supported, no further action needed
650  break;
651  } // switch(fAngleDist)
652 
653  delete histFile;
654  }
655 
656  //____________________________________________________________________________
657  bool SingleGen::PadVector(std::vector<double>& vec)
658  {
659  // check if the vec has the same size as fPDG
660  if (vec.size() != fPDG.size()) {
661  // if not padding out the vectors always cause an
662  // exception to be thrown if the vector in question
663  // is not the same size as the fPDG vector
664  // the exception is thrown in the reconfigure method
665  // that calls this one
666  if (!fPadOutVectors)
667  return false;
668  else if (fPadOutVectors) {
669  // if padding of vectors is desired but the vector in
670  // question has more than one entry it isn't clear
671  // what the padded values should be so cause
672  // an exception
673  if (vec.size() != 1) return false;
674 
675  // pad it out
676  vec.resize(fPDG.size(), vec[0]);
677 
678  } // end if padding out vectors
679  } // end if the vector size is not the same as fPDG
680 
681  return true;
682  }
683 
684  //____________________________________________________________________________
686  {
687 
689  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
690 
691  simb::MCTruth truth;
693  Sample(truth);
694 
695  MF_LOG_DEBUG("SingleGen") << truth;
696 
697  truthcol->push_back(truth);
698 
699  evt.put(std::move(truthcol));
700 
701  return;
702  }
703 
704  //____________________________________________________________________________
705  // Draw the type, momentum and position of a single particle from the
706  // FCIHL description
707  void SingleGen::SampleOne(unsigned int i, simb::MCTruth& mct)
708  {
709 
710  CLHEP::RandFlat flat(fEngine);
711  CLHEP::RandGaussQ gauss(fEngine);
712 
713  // Choose momentum
714  double p = 0.0;
715  double m = 0.0;
716  if (fPDist == kGAUS) { p = gauss.fire(fP0[i], fSigmaP[i]); }
717  else if (fPDist == kHIST) {
718  p = SelectFromHist(*(hPHist[i]));
719  }
720  else { // if (fPDist == kUNIF) {
721  p = fP0[i] + fSigmaP[i] * (2.0 * flat.fire() - 1.0);
722  }
723  // else {std::cout << "do not understand the value of PDist!";}
724 
725  static TDatabasePDG pdgt;
726  TParticlePDG* pdgp = pdgt.GetParticle(fPDG[i]);
727  if (pdgp) m = pdgp->Mass();
728 
729  // Choose position
730  TVector3 x;
731  if (fPosDist == kGAUS) {
732  x[0] = gauss.fire(fX0[i], fSigmaX[i]);
733  ;
734  x[1] = gauss.fire(fY0[i], fSigmaY[i]);
735  x[2] = gauss.fire(fZ0[i], fSigmaZ[i]);
736  }
737  else {
738  x[0] = fX0[i] + fSigmaX[i] * (2.0 * flat.fire() - 1.0);
739  x[1] = fY0[i] + fSigmaY[i] * (2.0 * flat.fire() - 1.0);
740  x[2] = fZ0[i] + fSigmaZ[i] * (2.0 * flat.fire() - 1.0);
741  }
742 
743  double t = 0.;
744  if (fTDist == kGAUS) { t = gauss.fire(fT0[i], fSigmaT[i]); }
745  else {
746  t = fT0[i] + fSigmaT[i] * (2.0 * flat.fire() - 1.0);
747  }
748 
749  TLorentzVector pos(x[0], x[1], x[2], t);
750 
751  // Choose angles
752  double thxz = 0;
753  double thyz = 0;
754 
755  double thyzrads = 0;
756  double thyzradsplussigma = 0;
757  double thyzradsminussigma = 0;
758 
759  if (fAngleDist == kGAUS) {
760  thxz = gauss.fire(fTheta0XZ[i], fSigmaThetaXZ[i]);
761  thyz = gauss.fire(fTheta0YZ[i], fSigmaThetaYZ[i]);
762  }
763  else if (fAngleDist == kHIST) { // Select thetaxz and thetayz from histogram
764  double thetaxz = 0;
765  double thetayz = 0;
766  SelectFromHist(*(hThetaXzYzHist[i]), thetaxz, thetayz);
767  thxz = (180. / M_PI) * thetaxz;
768  thyz = (180. / M_PI) * thetayz;
769  }
770  else {
771 
772  // Choose angles flat in phase space, which is flat in theta_xz
773  // and flat in sin(theta_yz).
774 
775  thxz = fTheta0XZ[i] + fSigmaThetaXZ[i] * (2.0 * flat.fire() - 1.0);
776 
777  thyzrads = std::asin(std::sin(
778  (M_PI / 180.) *
779  (fTheta0YZ
780  [i]))); //Taking asin of sin gives value between -Pi/2 and Pi/2 regardless of user input
781  thyzradsplussigma =
782  TMath::Min((thyzrads + ((M_PI / 180.) * fabs(fSigmaThetaYZ[i]))), M_PI / 2.);
783  thyzradsminussigma =
784  TMath::Max((thyzrads - ((M_PI / 180.) * fabs(fSigmaThetaYZ[i]))), -M_PI / 2.);
785 
786  //uncomment line to print angular variation info
787  //std::cout << "Central angle: " << (180./M_PI)*thyzrads << " Max angle: " << (180./M_PI)*thyzradsplussigma << " Min angle: " << (180./M_PI)*thyzradsminussigma << std::endl;
788 
789  double sinthyzmin = std::sin(thyzradsminussigma);
790  double sinthyzmax = std::sin(thyzradsplussigma);
791  double sinthyz = sinthyzmin + flat.fire() * (sinthyzmax - sinthyzmin);
792  thyz = (180. / M_PI) * std::asin(sinthyz);
793  }
794 
795  double thxzrad = thxz * M_PI / 180.0;
796  double thyzrad = thyz * M_PI / 180.0;
797 
798  TLorentzVector pvec(p * std::cos(thyzrad) * std::sin(thxzrad),
799  p * std::sin(thyzrad),
800  p * std::cos(thxzrad) * std::cos(thyzrad),
801  std::sqrt(p * p + m * m));
802 
803  // set track id to -i as these are all primary particles and have id <= 0
804  int trackid = -1 * (i + 1);
805  std::string primary("primary");
806 
807  simb::MCParticle part(trackid, fPDG[i], primary);
808  part.AddTrajectoryPoint(pos, pvec);
809 
810  //std::cout << "Px: " << pvec.Px() << " Py: " << pvec.Py() << " Pz: " << pvec.Pz() << std::endl;
811  //std::cout << "x: " << pos.X() << " y: " << pos.Y() << " z: " << pos.Z() << " time: " << pos.T() << std::endl;
812  //std::cout << "YZ Angle: " << (thyzrad * (180./M_PI)) << " XZ Angle: " << (thxzrad * (180./M_PI)) << std::endl;
813 
814  mct.Add(part);
815  }
816 
817  //____________________________________________________________________________
818  // Draw the type, momentum and position for all particles from the
819  // FCIHL description. Start positions will all match but momenta and angles drawn from
820  // distributions defined in the fhicls
822  {
823 
824  CLHEP::RandFlat flat(fEngine);
825  CLHEP::RandGaussQ gauss(fEngine);
826 
827  // Choose position
828  TVector3 x;
829  if (fPosDist == kGAUS) {
830  x[0] = gauss.fire(fX0[0], fSigmaX[0]);
831  ;
832  x[1] = gauss.fire(fY0[0], fSigmaY[0]);
833  x[2] = gauss.fire(fZ0[0], fSigmaZ[0]);
834  }
835  else {
836  x[0] = fX0[0] + fSigmaX[0] * (2.0 * flat.fire() - 1.0);
837  x[1] = fY0[0] + fSigmaY[0] * (2.0 * flat.fire() - 1.0);
838  x[2] = fZ0[0] + fSigmaZ[0] * (2.0 * flat.fire() - 1.0);
839  }
840 
841  double t = 0.;
842  if (fTDist == kGAUS) { t = gauss.fire(fT0[0], fSigmaT[0]); }
843  else {
844  t = fT0[0] + fSigmaT[0] * (2.0 * flat.fire() - 1.0);
845  }
846 
847  TLorentzVector pos(x[0], x[1], x[2], t);
848 
849  // loop through particles and select momenta and angles
850  for (unsigned int i(0); i < fPDG.size(); ++i) {
851  // Choose momentum
852  double p = 0.0;
853  double m = 0.0;
854  if (fPDist == kGAUS) { p = gauss.fire(fP0[i], fSigmaP[i]); }
855  else if (fPDist == kHIST) {
856  p = SelectFromHist(*(hPHist[i]));
857  }
858  else {
859  p = fP0[i] + fSigmaP[i] * (2.0 * flat.fire() - 1.0);
860  }
861 
862  static TDatabasePDG pdgt;
863  TParticlePDG* pdgp = pdgt.GetParticle(fPDG[i]);
864  if (pdgp) m = pdgp->Mass();
865 
866  // Choose angles
867  double thxz = 0;
868  double thyz = 0;
869 
870  double thyzrads = 0;
871  double thyzradsplussigma = 0;
872  double thyzradsminussigma = 0;
873 
874  if (fAngleDist == kGAUS) {
875  thxz = gauss.fire(fTheta0XZ[i], fSigmaThetaXZ[i]);
876  thyz = gauss.fire(fTheta0YZ[i], fSigmaThetaYZ[i]);
877  }
878  else if (fAngleDist == kHIST) {
879  double thetaxz = 0;
880  double thetayz = 0;
881  SelectFromHist(*(hThetaXzYzHist[i]), thetaxz, thetayz);
882  thxz = (180. / M_PI) * thetaxz;
883  thyz = (180. / M_PI) * thetayz;
884  }
885  else {
886 
887  // Choose angles flat in phase space, which is flat in theta_xz
888  // and flat in sin(theta_yz).
889 
890  thxz = fTheta0XZ[i] + fSigmaThetaXZ[i] * (2.0 * flat.fire() - 1.0);
891 
892  thyzrads = std::asin(std::sin(
893  (M_PI / 180.) *
894  (fTheta0YZ
895  [i]))); //Taking asin of sin gives value between -Pi/2 and Pi/2 regardless of user input
896  thyzradsplussigma =
897  TMath::Min((thyzrads + ((M_PI / 180.) * fabs(fSigmaThetaYZ[i]))), M_PI / 2.);
898  thyzradsminussigma =
899  TMath::Max((thyzrads - ((M_PI / 180.) * fabs(fSigmaThetaYZ[i]))), -M_PI / 2.);
900 
901  //uncomment line to print angular variation info
902  //std::cout << "Central angle: " << (180./M_PI)*thyzrads << " Max angle: " << (180./M_PI)*thyzradsplussigma << " Min angle: " << (180./M_PI)*thyzradsminussigma << std::endl;
903 
904  double sinthyzmin = std::sin(thyzradsminussigma);
905  double sinthyzmax = std::sin(thyzradsplussigma);
906  double sinthyz = sinthyzmin + flat.fire() * (sinthyzmax - sinthyzmin);
907  thyz = (180. / M_PI) * std::asin(sinthyz);
908  }
909 
910  double thxzrad = thxz * M_PI / 180.0;
911  double thyzrad = thyz * M_PI / 180.0;
912 
913  TLorentzVector pvec(p * std::cos(thyzrad) * std::sin(thxzrad),
914  p * std::sin(thyzrad),
915  p * std::cos(thxzrad) * std::cos(thyzrad),
916  std::sqrt(p * p + m * m));
917 
918  // set track id to -i as these are all primary particles and have id <= 0
919  int trackid = -1 * (i + 1);
920  std::string primary("primary");
921 
922  simb::MCParticle part(trackid, fPDG[i], primary);
923  part.AddTrajectoryPoint(pos, pvec);
924 
925  //std::cout << "Px: " << pvec.Px() << " Py: " << pvec.Py() << " Pz: " << pvec.Pz() << std::endl;
926  //std::cout << "x: " << pos.X() << " y: " << pos.Y() << " z: " << pos.Z() << " time: " << pos.T() << std::endl;
927  //std::cout << "YZ Angle: " << (thyzrad * (180./M_PI)) << " XZ Angle: " << (thxzrad * (180./M_PI)) << std::endl;
928  mct.Add(part);
929  }
930  }
931 
932  //____________________________________________________________________________
934  {
935 
936  switch (fMode) {
937  case 0: // List generation mode: every event will have one of each
938  // particle species in the fPDG array
939  if (fSingleVertex) { SampleMany(mct); }
940  else {
941  for (unsigned int i = 0; i < fPDG.size(); ++i) {
942  SampleOne(i, mct);
943  } //end loop over particles
944  }
945  break;
946  case 1: // Random selection mode: every event will exactly one particle
947  // selected randomly from the fPDG array
948  {
949  CLHEP::RandFlat flat(fEngine);
950 
951  unsigned int i = flat.fireInt(fPDG.size());
952  SampleOne(i, mct);
953  } break;
954  default:
955  mf::LogWarning("UnrecognizeOption")
956  << "SingleGen does not recognize ParticleSelectionMode " << fMode;
957  break;
958  } // switch on fMode
959 
960  return;
961  }
962 
963  //____________________________________________________________________________
964  void SingleGen::printVecs(std::vector<std::string> const& list)
965  {
966 
967  mf::LogInfo("SingleGen") << " You are using vector values for SingleGen configuration.\n "
968  << " Some of the configuration vectors may have been padded out ,"
969  << " because they (weren't) as long as the pdg vector"
970  << " in your configuration. \n"
971  << " The new input particle configuration is:\n";
972 
973  std::string values;
974  for (size_t i = 0; i <= 1; ++i) { // list.size(); ++i){
975 
976  values.append(list[i]);
977  values.append(": [ ");
978 
979  for (size_t e = 0; e < fPDG.size(); ++e) {
980  std::stringstream buf;
981  buf.width(10);
982  if (i == 0) buf << fPDG[e] << ", ";
983  buf.precision(5);
984  if (i == 1) buf << fP0[e] << ", ";
985  if (i == 2) buf << fSigmaP[e] << ", ";
986  if (i == 3) buf << fX0[e] << ", ";
987  if (i == 4) buf << fY0[e] << ", ";
988  if (i == 5) buf << fZ0[e] << ", ";
989  if (i == 6) buf << fSigmaX[e] << ", ";
990  if (i == 7) buf << fSigmaY[e] << ", ";
991  if (i == 8) buf << fSigmaZ[e] << ", ";
992  if (i == 9) buf << fTheta0XZ[e] << ", ";
993  if (i == 10) buf << fTheta0YZ[e] << ", ";
994  if (i == 11) buf << fSigmaThetaXZ[e] << ", ";
995  if (i == 12) buf << fSigmaThetaYZ[e] << ", ";
996  if (i == 13) buf << fT0[e] << ", ";
997  if (i == 14) buf << fSigmaT[e] << ", ";
998  values.append(buf.str());
999  }
1000 
1001  values.erase(values.find_last_of(","));
1002  values.append(" ] \n");
1003 
1004  } // end loop over vector names in list
1005 
1006  mf::LogInfo("SingleGen") << values;
1007 
1008  return;
1009  }
1010 
1011  //____________________________________________________________________________
1012  double SingleGen::SelectFromHist(const TH1& h) // select from a 1D histogram
1013  {
1014  CLHEP::RandFlat flat(fEngine);
1015 
1016  double throw_value = h.Integral() * flat.fire();
1017  double cum_value(0);
1018  for (int i(0); i < h.GetNbinsX() + 1; ++i) {
1019  cum_value += h.GetBinContent(i);
1020  if (throw_value < cum_value) { return flat.fire() * h.GetBinWidth(i) + h.GetBinLowEdge(i); }
1021  }
1022  return throw_value; // for some reason we've gone through all bins and failed?
1023  }
1024  //____________________________________________________________________________
1025  void SingleGen::SelectFromHist(const TH2& h, double& x, double& y) // select from a 2D histogram
1026  {
1027  CLHEP::RandFlat flat(fEngine);
1028 
1029  double throw_value = h.Integral() * flat.fire();
1030  double cum_value(0);
1031  for (int i(0); i < h.GetNbinsX() + 1; ++i) {
1032  for (int j(0); j < h.GetNbinsY() + 1; ++j) {
1033  cum_value += h.GetBinContent(i, j);
1034  if (throw_value < cum_value) {
1035  x = flat.fire() * h.GetXaxis()->GetBinWidth(i) + h.GetXaxis()->GetBinLowEdge(i);
1036  y = flat.fire() * h.GetYaxis()->GetBinWidth(j) + h.GetYaxis()->GetBinLowEdge(j);
1037  return;
1038  }
1039  }
1040  }
1041  return; // for some reason we've gone through all bins and failed?
1042  }
1043  //____________________________________________________________________________
1044 
1045 } //end namespace evgen
1046 
1047 namespace evgen {
1048 
1050 
1051 } //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.
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:258
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 (s)
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
fhicl::OptionalAtom< rndm::NuRandomService::seed_t > Seed
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 detector geometry.
fhicl::Sequence< double > T0
fhicl::Atom< std::string > HistogramFile
std::vector< double > fT0
Central t position (s) 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
long seed
Definition: chem4.cc:67
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
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 geometry of one entire detector.
Definition: GeometryCore.h:119
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
#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:203
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
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 .
art::detail::EngineCreator::seed_t seed_t
std::vector< double > fSigmaY
Variation in y position (cm)