13 #include <initializer_list> 26 #include "cetlib/filesystem.h" 27 #include "cetlib/search_path.h" 28 #include "cetlib_except/exception.h" 49 #include "TDatabasePDG.h" 55 #include "CLHEP/Random/RandFlat.h" 56 #include "CLHEP/Random/RandGaussQ.h" 69 Name(
"ParticleSelectionMode"),
70 Comment(
"generate one particle, or one particle per PDG ID: " +
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")};
79 Comment(
"PDG ID of the particles to be generated; this is the key " 80 "for the other options marked as \"per PDG ID\"")};
88 Comment(
"central momentum (GeV/c) to generate"),
92 Comment(
"variation in momenta (GeV/c)"),
97 Comment(
"central x position (cm) in world coordinates [per PDG ID]")};
101 Comment(
"central y position (cm) in world coordinates [per PDG ID]")};
105 Comment(
"central z position (cm) in world coordinates [per PDG ID]")};
111 Comment(
"variation (radius or RMS) in x position (cm) [per PDG ID]")};
115 Comment(
"variation (radius or RMS) in y position (cm) [per PDG ID]")};
119 Comment(
"variation (radius or RMS) in z position (cm) [per PDG ID]")};
123 Comment(
"variation (semi-interval or RMS) in time (ns) [per PDG ID]")};
126 Comment(
"distribution of starting position: " +
134 Name(
"SingleVertex"),
135 Comment(
"if true, all particles are produced at the same location"),
139 Comment(
"angle from Z axis on world X-Z plane (degrees)")};
142 Comment(
"angle from Z axis on world Y-Z plane (degrees)")};
145 Comment(
"variation in angle in X-Z plane (degrees)")};
148 Comment(
"variation in angle in Y-Z plane (degrees)")};
156 Name(
"HistogramFile"),
157 Comment(
"ROOT file containing the required distributions for the generation"),
162 Comment(
"name of the histograms of momentum distributions"),
173 Name(
"ThetaXzYzHist"),
174 Comment(
"name of the histograms of angular (X-Z and Y-Z) distribution"),
179 Comment(
"override the random number generator seed"),
204 void printVecs(std::vector<std::string>
const& list);
205 bool PadVector(std::vector<double>& vec);
233 std::vector<int>
fPDG;
260 std::vector<std::unique_ptr<TH1>>
hPHist;
261 std::vector<std::unique_ptr<TH2>>
303 template <
typename OptionList>
304 static auto selectOption(std::string Option, OptionList
const& allowedOptions)
320 template <
typename OptionList>
322 OptionList
const& allowedOptions,
324 std::initializer_list<typename OptionList::value_type::first_type> exclude);
326 template <
typename OptionList>
327 static std::string
presentOptions(OptionList
const& allowedOptions,
bool printKey =
true)
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>");
345 std::map<int, std::string> names;
353 std::map<int, std::string> names;
354 names[int(
kUNIF)] =
"uniform";
355 names[int(
kGAUS)] =
"Gaussian";
356 names[int(
kHIST)] =
"histograms";
365 template <
typename OptionList>
369 using key_type =
typename OptionList::value_type::first_type;
370 using tolower_type = int (*)(int);
371 auto toLower = [](
auto const& S) {
374 std::transform(S.cbegin(), S.cend(), std::back_inserter(s), (tolower_type)&std::tolower);
377 auto option = toLower(Option);
378 for (
auto const& candidate : allowedOptions) {
379 if (toLower(candidate.second) == option)
return candidate.first;
383 key_type num = std::stoi(Option, &end);
384 if (allowedOptions.count(num) && (end == Option.length()))
return num;
386 catch (std::invalid_argument
const&) {
388 throw std::runtime_error(
"Option '" + Option +
"' not supported.");
391 template <
typename OptionList>
393 OptionList
const& allowedOptions,
395 std::initializer_list<typename OptionList::value_type::first_type> exclude
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) +
'\"';
411 template <
typename OptionList>
413 OptionList
const& allowedOptions,
417 auto iOption = allowedOptions.find(optionKey);
418 return (iOption != allowedOptions.end()) ? iOption->second : defName;
429 constexpr
double kAlphaMass = 3.727379240;
430 if (particle_id == 1000020040) {
434 static TDatabasePDG pdgt;
435 TParticlePDG* pdgp = pdgt.GetParticle(particle_id);
436 return pdgp ? pdgp->Mass() : 0.;
441 if (particle_id == 1000020040) {
452 ,
fPDG(config().PDG())
473 ,
fPHist(config().PHist())
476 ->registerAndSeedEngine(
createEngine(config().Seed()),
"",
"", config().Seed()))
481 produces<std::vector<simb::MCTruth>>();
482 produces<sumdata::RunData, art::InRun>();
497 std::vector<std::string> vlist(15);
507 vlist[9] =
"Theta0XZ";
508 vlist[10] =
"Theta0YZ";
509 vlist[11] =
"SigmaThetaXZ";
510 vlist[12] =
"SigmaThetaYZ";
512 vlist[14] =
"SigmaT";
521 if (!this->
PadVector(
fP0)) { list.append(vlist[1].append(
", \n")); }
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")); }
534 if (!this->
PadVector(
fT0)) { list.append(vlist[13].append(
", \n")); }
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";
548 TFile* histFile =
nullptr;
554 if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
558 <<
"Cannot open ROOT file specified in parameter HistogramFile: \"" <<
fHistFileName 564 <<
"ROOT file specified in parameter HistogramFile: \"" <<
fHistFileName 565 <<
"\" does not exist!";
570 std::string relative_filename{
"./"};
572 if (cet::file_exists(relative_filename)) {
573 histFile =
new TFile(relative_filename.c_str());
574 if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
578 <<
"Cannot open ROOT file found using relative path and originally specified in " 579 "parameter HistogramFile: \"" 580 << relative_filename <<
'"';
584 cet::search_path sp{
"FW_SEARCH_PATH"};
585 std::string found_filename;
589 <<
"Cannot find ROOT file in current directory nor on FW_SEARCH_PATH specified in " 590 "parameter HistogramFile: \"" 593 histFile =
new TFile(found_filename.c_str());
594 if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
598 <<
"Cannot open ROOT file found on FW_SEARCH_PATH and originally specified in " 599 "parameter HistogramFile: \"" 600 << found_filename <<
'"';
637 <<
fPHist.size() <<
" momentum histograms to describe " <<
fPDG.size()
638 <<
" particle types...";
641 for (
auto const& histName :
fPHist) {
642 TH1* pHist =
dynamic_cast<TH1*
>(histFile->Get(histName.c_str()));
645 <<
"Failed to read momentum histogram '" << histName <<
"' from '" 646 << histFile->GetPath() <<
"\'";
648 pHist->SetDirectory(
nullptr);
649 hPHist.emplace_back(pHist);
661 <<
" particle types...";
665 TH2* pHist =
dynamic_cast<TH2*
>(histFile->Get(histName.c_str()));
668 <<
"Failed to read direction histogram '" << histName <<
"' from '" 669 << histFile->GetPath() <<
"\'";
671 pHist->SetDirectory(
nullptr);
685 if (vec.size() !=
fPDG.size()) {
698 if (vec.size() != 1)
return false;
701 vec.resize(
fPDG.size(), vec[0]);
714 std::unique_ptr<std::vector<simb::MCTruth>> truthcol(
new std::vector<simb::MCTruth>);
722 truthcol->push_back(truth);
724 evt.
put(std::move(truthcol));
736 CLHEP::RandGaussQ gauss(
fEngine);
746 p =
fP0[i] +
fSigmaP[i] * (2.0 * flat.fire() - 1.0);
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);
768 t =
fT0[i] +
fSigmaT[i] * (2.0 * flat.fire() - 1.0);
771 TLorentzVector pos(x[0], x[1], x[2], t);
778 double thyzradsplussigma = 0;
779 double thyzradsminussigma = 0;
789 thxz = (180. / M_PI) * thetaxz;
790 thyz = (180. / M_PI) * thetayz;
799 thyzrads = std::asin(std::sin(
804 TMath::Min((thyzrads + ((M_PI / 180.) * fabs(
fSigmaThetaYZ[i]))), M_PI / 2.);
806 TMath::Max((thyzrads - ((M_PI / 180.) * fabs(
fSigmaThetaYZ[i]))), -M_PI / 2.);
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);
817 double thxzrad = thxz * M_PI / 180.0;
818 double thyzrad = thyz * M_PI / 180.0;
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));
826 int trackid = -1 * (i + 1);
827 std::string primary(
"primary");
830 part.AddTrajectoryPoint(pos, pvec);
846 CLHEP::RandGaussQ gauss(
fEngine);
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);
865 t =
fT0[0] +
fSigmaT[0] * (2.0 * flat.fire() - 1.0);
868 TLorentzVector pos(x[0], x[1], x[2], t);
871 for (
unsigned int i(0); i <
fPDG.size(); ++i) {
880 p =
fP0[i] +
fSigmaP[i] * (2.0 * flat.fire() - 1.0);
890 double thyzradsplussigma = 0;
891 double thyzradsminussigma = 0;
901 thxz = (180. / M_PI) * thetaxz;
902 thyz = (180. / M_PI) * thetayz;
911 thyzrads = std::asin(std::sin(
916 TMath::Min((thyzrads + ((M_PI / 180.) * fabs(
fSigmaThetaYZ[i]))), M_PI / 2.);
918 TMath::Max((thyzrads - ((M_PI / 180.) * fabs(
fSigmaThetaYZ[i]))), -M_PI / 2.);
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);
929 double thxzrad = thxz * M_PI / 180.0;
930 double thyzrad = thyz * M_PI / 180.0;
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));
938 int trackid = -1 * (i + 1);
939 std::string primary(
"primary");
942 part.AddTrajectoryPoint(pos, pvec);
960 for (
unsigned int i = 0; i <
fPDG.size(); ++i) {
970 unsigned int i = flat.fireInt(
fPDG.size());
975 <<
"SingleGen does not recognize ParticleSelectionMode " <<
fMode;
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";
993 for (
size_t i = 0; i <= 1; ++i) {
995 values.append(list[i]);
996 values.append(
": [ ");
998 for (
size_t e = 0;
e <
fPDG.size(); ++
e) {
999 std::stringstream buf;
1001 if (i == 0) buf <<
fPDG[
e] <<
", ";
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] <<
", ";
1015 if (i == 13) buf <<
fT0[
e] <<
", ";
1016 if (i == 14) buf <<
fSigmaT[
e] <<
", ";
1017 values.append(buf.str());
1020 values.erase(values.find_last_of(
","));
1021 values.append(
" ] \n");
1033 CLHEP::RandFlat flat(
fEngine);
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); }
1046 CLHEP::RandFlat flat(
fEngine);
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);
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.
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)
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
EDProducer(fhicl::ParameterSet const &pset)
std::vector< std::unique_ptr< TH1 > > hPHist
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
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
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
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={})
static constexpr int kSelectOneRandPart
fhicl::Sequence< double > Theta0YZ
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
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.
fhicl::Atom< std::string > AngleDist
void produce(art::Event &evt) override
#define DEFINE_ART_MODULE(klass)
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
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.
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
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)
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.
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.
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.
Event generator information.
void beginRun(art::Run &run) override
Act on begin of run: write "RunData" information (sumdata::RunData).
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
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)