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 (s) [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")};
203 void printVecs(std::vector<std::string>
const& list);
204 bool PadVector(std::vector<double>& vec);
228 std::vector<int>
fPDG;
255 std::vector<std::unique_ptr<TH1>>
hPHist;
256 std::vector<std::unique_ptr<TH2>>
298 template <
typename OptionList>
299 static auto selectOption(std::string Option, OptionList
const& allowedOptions)
315 template <
typename OptionList>
317 OptionList
const& allowedOptions,
319 std::initializer_list<typename OptionList::value_type::first_type> exclude);
321 template <
typename OptionList>
322 static std::string
presentOptions(OptionList
const& allowedOptions,
bool printKey =
true)
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>");
340 std::map<int, std::string> names;
348 std::map<int, std::string> names;
349 names[int(
kUNIF)] =
"uniform";
350 names[int(
kGAUS)] =
"Gaussian";
351 names[int(
kHIST)] =
"histograms";
360 template <
typename OptionList>
364 using key_type =
typename OptionList::value_type::first_type;
365 using tolower_type = int (*)(int);
366 auto toLower = [](
auto const& S) {
369 std::transform(S.cbegin(), S.cend(), std::back_inserter(s), (tolower_type)&std::tolower);
372 auto option = toLower(Option);
373 for (
auto const& candidate : allowedOptions) {
374 if (toLower(candidate.second) == option)
return candidate.first;
378 key_type num = std::stoi(Option, &end);
379 if (allowedOptions.count(num) && (end == Option.length()))
return num;
381 catch (std::invalid_argument
const&) {
383 throw std::runtime_error(
"Option '" + Option +
"' not supported.");
386 template <
typename OptionList>
388 OptionList
const& allowedOptions,
390 std::initializer_list<typename OptionList::value_type::first_type> exclude
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) +
'\"';
406 template <
typename OptionList>
408 OptionList
const& allowedOptions,
412 auto iOption = allowedOptions.find(optionKey);
413 return (iOption != allowedOptions.end()) ? iOption->second : defName;
427 ,
fPDG(config().PDG())
448 ,
fPHist(config().PHist())
454 if (config().Seed(seed)) {
fEngine.setSeed(seed, 0 ); }
456 produces<std::vector<simb::MCTruth>>();
457 produces<sumdata::RunData, art::InRun>();
472 std::vector<std::string> vlist(15);
482 vlist[9] =
"Theta0XZ";
483 vlist[10] =
"Theta0YZ";
484 vlist[11] =
"SigmaThetaXZ";
485 vlist[12] =
"SigmaThetaYZ";
487 vlist[14] =
"SigmaT";
496 if (!this->
PadVector(
fP0)) { list.append(vlist[1].append(
", \n")); }
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")); }
509 if (!this->
PadVector(
fT0)) { list.append(vlist[13].append(
", \n")); }
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";
523 TFile* histFile =
nullptr;
529 if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
533 <<
"Cannot open ROOT file specified in parameter HistogramFile: \"" <<
fHistFileName 539 <<
"ROOT file specified in parameter HistogramFile: \"" <<
fHistFileName 540 <<
"\" does not exist!";
545 std::string relative_filename{
"./"};
547 if (cet::file_exists(relative_filename)) {
548 histFile =
new TFile(relative_filename.c_str());
549 if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
553 <<
"Cannot open ROOT file found using relative path and originally specified in " 554 "parameter HistogramFile: \"" 555 << relative_filename <<
'"';
559 cet::search_path sp{
"FW_SEARCH_PATH"};
560 std::string found_filename;
564 <<
"Cannot find ROOT file in current directory nor on FW_SEARCH_PATH specified in " 565 "parameter HistogramFile: \"" 568 histFile =
new TFile(found_filename.c_str());
569 if (!histFile || histFile->IsZombie() || !histFile->IsOpen()) {
573 <<
"Cannot open ROOT file found on FW_SEARCH_PATH and originally specified in " 574 "parameter HistogramFile: \"" 575 << found_filename <<
'"';
612 <<
fPHist.size() <<
" momentum histograms to describe " <<
fPDG.size()
613 <<
" particle types...";
616 for (
auto const& histName :
fPHist) {
617 TH1* pHist =
dynamic_cast<TH1*
>(histFile->Get(histName.c_str()));
620 <<
"Failed to read momentum histogram '" << histName <<
"' from '" 621 << histFile->GetPath() <<
"\'";
623 pHist->SetDirectory(
nullptr);
624 hPHist.emplace_back(pHist);
636 <<
" particle types...";
640 TH2* pHist =
dynamic_cast<TH2*
>(histFile->Get(histName.c_str()));
643 <<
"Failed to read direction histogram '" << histName <<
"' from '" 644 << histFile->GetPath() <<
"\'";
646 pHist->SetDirectory(
nullptr);
660 if (vec.size() !=
fPDG.size()) {
673 if (vec.size() != 1)
return false;
676 vec.resize(
fPDG.size(), vec[0]);
689 std::unique_ptr<std::vector<simb::MCTruth>> truthcol(
new std::vector<simb::MCTruth>);
697 truthcol->push_back(truth);
699 evt.
put(std::move(truthcol));
711 CLHEP::RandGaussQ gauss(
fEngine);
721 p =
fP0[i] +
fSigmaP[i] * (2.0 * flat.fire() - 1.0);
725 static TDatabasePDG pdgt;
726 TParticlePDG* pdgp = pdgt.GetParticle(
fPDG[i]);
727 if (pdgp) m = pdgp->Mass();
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);
746 t =
fT0[i] +
fSigmaT[i] * (2.0 * flat.fire() - 1.0);
749 TLorentzVector pos(x[0], x[1], x[2], t);
756 double thyzradsplussigma = 0;
757 double thyzradsminussigma = 0;
767 thxz = (180. / M_PI) * thetaxz;
768 thyz = (180. / M_PI) * thetayz;
777 thyzrads = std::asin(std::sin(
782 TMath::Min((thyzrads + ((M_PI / 180.) * fabs(
fSigmaThetaYZ[i]))), M_PI / 2.);
784 TMath::Max((thyzrads - ((M_PI / 180.) * fabs(
fSigmaThetaYZ[i]))), -M_PI / 2.);
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);
795 double thxzrad = thxz * M_PI / 180.0;
796 double thyzrad = thyz * M_PI / 180.0;
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));
804 int trackid = -1 * (i + 1);
805 std::string primary(
"primary");
825 CLHEP::RandGaussQ gauss(
fEngine);
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);
844 t =
fT0[0] +
fSigmaT[0] * (2.0 * flat.fire() - 1.0);
847 TLorentzVector pos(x[0], x[1], x[2], t);
850 for (
unsigned int i(0); i <
fPDG.size(); ++i) {
859 p =
fP0[i] +
fSigmaP[i] * (2.0 * flat.fire() - 1.0);
862 static TDatabasePDG pdgt;
863 TParticlePDG* pdgp = pdgt.GetParticle(
fPDG[i]);
864 if (pdgp) m = pdgp->Mass();
871 double thyzradsplussigma = 0;
872 double thyzradsminussigma = 0;
882 thxz = (180. / M_PI) * thetaxz;
883 thyz = (180. / M_PI) * thetayz;
892 thyzrads = std::asin(std::sin(
897 TMath::Min((thyzrads + ((M_PI / 180.) * fabs(
fSigmaThetaYZ[i]))), M_PI / 2.);
899 TMath::Max((thyzrads - ((M_PI / 180.) * fabs(
fSigmaThetaYZ[i]))), -M_PI / 2.);
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);
910 double thxzrad = thxz * M_PI / 180.0;
911 double thyzrad = thyz * M_PI / 180.0;
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));
919 int trackid = -1 * (i + 1);
920 std::string primary(
"primary");
941 for (
unsigned int i = 0; i <
fPDG.size(); ++i) {
951 unsigned int i = flat.fireInt(
fPDG.size());
956 <<
"SingleGen does not recognize ParticleSelectionMode " <<
fMode;
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";
974 for (
size_t i = 0; i <= 1; ++i) {
976 values.append(list[i]);
977 values.append(
": [ ");
979 for (
size_t e = 0;
e <
fPDG.size(); ++
e) {
980 std::stringstream buf;
982 if (i == 0) buf <<
fPDG[
e] <<
", ";
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] <<
", ";
996 if (i == 13) buf <<
fT0[
e] <<
", ";
997 if (i == 14) buf <<
fSigmaT[
e] <<
", ";
998 values.append(buf.str());
1001 values.erase(values.find_last_of(
","));
1002 values.append(
" ] \n");
1014 CLHEP::RandFlat flat(
fEngine);
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); }
1027 CLHEP::RandFlat flat(
fEngine);
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);
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)
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 (s)
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={})
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
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 detector geometry.
fhicl::Sequence< double > T0
fhicl::Atom< std::string > HistogramFile
std::vector< double > fT0
Central t position (s) 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
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.
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
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
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 .
art::detail::EngineCreator::seed_t seed_t
std::vector< double > fSigmaY
Variation in y position (cm)