53 #include "cetlib/exempt_ptr.h" 54 #include "cetlib/search_path.h" 55 #include "cetlib_except/exception.h" 89 #include "TGenPhaseSpace.h" 91 #include "TGeoManager.h" 92 #include "TGeoMaterial.h" 94 #include "TGeoVolume.h" 98 #include "TLorentzVector.h" 101 #include "CLHEP/Random/RandFlat.h" 102 #include "CLHEP/Random/RandPoisson.h" 209 std::size_t addvolume(std::string
const& volumeName);
214 TLorentzVector dirCalc(
double p,
double m);
216 void readfile(std::string nuclide, std::string
const& filename, std::string
const& PathToFile);
217 void samplespectrum(std::string nuclide,
int& itype,
double& t,
double& m,
double& p);
219 void Ar42Gamma2(
std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>& v_prods);
220 void Ar42Gamma3(
std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>& v_prods);
221 void Ar42Gamma4(
std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>& v_prods);
222 void Ar42Gamma5(
std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>& v_prods);
225 double samplefromth1d(TH1D&
hist);
228 template <
typename Stream>
229 void dumpNuclideSettings(Stream&& out,
231 std::string
const& volumeName = {})
const;
234 std::pair<double, double> defaulttimewindow()
const;
247 std::vector<std::string>
280 constexpr
double m_e = 0.000510998928;
281 constexpr
double m_alpha = 3.727379240;
282 constexpr
double m_neutron = 0.9395654133;
290 ,
fX0{pset.get<std::vector<double>>(
"X0", {})}
291 ,
fY0{pset.get<std::vector<double>>(
"Y0", {})}
292 ,
fZ0{pset.get<std::vector<double>>(
"Z0", {})}
293 ,
fX1{pset.get<std::vector<double>>(
"X1", {})}
294 ,
fY1{pset.get<std::vector<double>>(
"Y1", {})}
295 ,
fZ1{pset.get<std::vector<double>>(
"Z1", {})}
303 produces<std::vector<simb::MCTruth>>();
304 produces<sumdata::RunData, art::InRun>();
306 auto const fPathToFile = pset.get<std::string>(
"PathToFile",
"Radionuclides/");
307 auto const nuclide = pset.get<std::vector<std::string>>(
"Nuclide");
308 auto const material = pset.get<std::vector<std::string>>(
"Material");
309 auto const Bq = pset.get<std::vector<double>>(
"BqPercc");
310 auto t0 = pset.get<std::vector<double>>(
"T0", {});
311 auto t1 = pset.get<std::vector<double>>(
"T1", {});
313 if (fPathToFile.find_last_of(
'/') !=
314 fPathToFile.length() - 1) {
316 <<
"Last entry of the path to radiological files needs to be a forward slash ('/')\n";
319 if (
fT0.empty() ||
fT1.empty()) {
320 if (!
fT0.empty() || !
fT1.empty()) {
322 <<
"RadioGen T0 and T1 need to be both non-empty, or both empty" 324 <<
fT0.size() <<
" entries and T1 has " <<
fT0.size() <<
")\n";
327 t0.push_back(defaultT0);
328 t1.push_back(defaultT1);
337 mf::LogInfo(
"RadioGen") <<
"Configuring activity:";
340 fNuclide.push_back(nuclide.at(iCoord));
341 fMaterial.push_back(material.at(iCoord));
342 fBq.push_back(Bq.at(iCoord));
344 fT0.push_back(t0[std::min(iCoord, t0.size() - 1U)]);
345 fT1.push_back(t1[std::min(iCoord, t1.size() - 1U)]);
351 std::size_t iNext =
fX0.size();
352 auto const volumes = pset.get<std::vector<std::string>>(
"Volumes", {});
355 auto const nVolumes =
addvolume(volName);
358 <<
"No volume named '" << volName <<
"' was found!\n";
361 std::size_t
const iVolFirst =
fNuclide.size();
363 fNuclide.push_back(nuclide.at(iNext));
365 fBq.push_back(Bq.at(iNext));
367 fT0.push_back(t0[std::min(iNext, t0.size() - 1U)]);
368 fT1.push_back(t1[std::min(iNext, t1.size() - 1U)]);
376 unsigned int nsize =
fNuclide.size();
382 throw cet::exception(
"RadioGen") <<
"Different size Material vector and Nuclide vector\n";
383 if (
fBq.size() != nsize)
384 throw cet::exception(
"RadioGen") <<
"Different size Bq vector and Nuclide vector\n";
385 if (
fT0.size() != nsize)
386 throw cet::exception(
"RadioGen") <<
"Different size T0 vector and Nuclide vector\n";
387 if (
fT1.size() != nsize)
388 throw cet::exception(
"RadioGen") <<
"Different size T1 vector and Nuclide vector\n";
389 if (
fX0.size() != nsize)
390 throw cet::exception(
"RadioGen") <<
"Different size X0 vector and Nuclide vector\n";
391 if (
fY0.size() != nsize)
392 throw cet::exception(
"RadioGen") <<
"Different size Y0 vector and Nuclide vector\n";
393 if (
fZ0.size() != nsize)
394 throw cet::exception(
"RadioGen") <<
"Different size Z0 vector and Nuclide vector\n";
395 if (
fX1.size() != nsize)
396 throw cet::exception(
"RadioGen") <<
"Different size X1 vector and Nuclide vector\n";
397 if (
fY1.size() != nsize)
398 throw cet::exception(
"RadioGen") <<
"Different size Y1 vector and Nuclide vector\n";
399 if (
fZ1.size() != nsize)
400 throw cet::exception(
"RadioGen") <<
"Different size Z1 vector and Nuclide vector\n";
402 for (std::string& nuclideName :
fNuclide) {
403 if (nuclideName ==
"39Ar") {
readfile(
"39Ar",
"Argon_39.root", fPathToFile); }
404 else if (nuclideName ==
"60Co") {
405 readfile(
"60Co",
"Cobalt_60.root", fPathToFile);
407 else if (nuclideName ==
"85Kr") {
408 readfile(
"85Kr",
"Krypton_85.root", fPathToFile);
410 else if (nuclideName ==
"40K") {
411 readfile(
"40K",
"Potassium_40.root", fPathToFile);
413 else if (nuclideName ==
"232Th") {
414 readfile(
"232Th",
"Thorium_232.root", fPathToFile);
416 else if (nuclideName ==
"238U") {
417 readfile(
"238U",
"Uranium_238.root", fPathToFile);
419 else if (nuclideName ==
"222Rn") {
422 else if (nuclideName ==
"59Ni") {
425 else if (nuclideName ==
"42Ar") {
442 readfile(
"42Ar_5",
"Argon_42_5.root", fPathToFile);
446 std::string searchName = nuclideName;
447 searchName +=
".root";
448 readfile(nuclideName, searchName, fPathToFile);
464 std::unique_ptr<std::vector<simb::MCTruth>> truthcol(
new std::vector<simb::MCTruth>);
470 for (
unsigned int i = 0; i <
fNuclide.size(); ++i) {
475 truthcol->push_back(truth);
476 evt.
put(std::move(truthcol));
483 auto const& geom = *(lar::providerFrom<geo::Geometry>());
485 std::vector<geo::GeoNodePath> volumePaths;
486 auto findVolume = [&volumePaths, volumeName](
auto& path) {
487 if (path.current().GetVolume()->GetName() == volumeName) volumePaths.push_back(path);
493 navigator.
apply(findVolume);
500 TGeoShape
const* pShape = path.current().GetVolume()->GetShape();
501 auto pBox =
dynamic_cast<TGeoBBox const*
>(pShape);
503 throw cet::exception(
"RadioGen") <<
"Volume '" << path.current().GetName() <<
"' is a " 504 << pShape->IsA()->GetName() <<
", not a TGeoBBox.\n";
507 geo::Point_t const origin = geo::vect::makeFromCoords<geo::Point_t>(pBox->GetOrigin());
508 geo::Vector_t const diag = {pBox->GetDX(), pBox->GetDY(), pBox->GetDZ()};
517 trans.Transform(origin - diag, min);
518 trans.Transform(origin + diag, max);
523 fX0.push_back(std::min(min.X(), max.X()));
524 fY0.push_back(std::min(min.Y(), max.Y()));
525 fZ0.push_back(std::min(min.Z(), max.Z()));
526 fX1.push_back(std::max(min.X(), max.X()));
527 fY1.push_back(std::max(min.Y(), max.Y()));
528 fZ1.push_back(std::max(min.Z(), max.Z()));
532 return volumePaths.size();
544 CLHEP::RandPoisson poisson(
fEngine);
552 long ndecays = poisson.shoot(rate);
554 std::regex
const re_material{
fMaterial[i]};
555 for (
unsigned int idecay = 0; idecay < ndecays; idecay++) {
562 fX0[i] + flat.fire() * (
fX1[i] -
fX0[i]),
563 fY0[i] + flat.fire() * (
fY1[i] -
fY0[i]),
564 fZ0[i] + flat.fire() * (
fZ1[i] -
fZ0[i]),
568 std::string volmaterial =
569 geomanager->FindNode(pos.X(), pos.Y(), pos.Z())->GetMedium()->GetMaterial()->GetName();
570 if (!std::regex_match(volmaterial, re_material))
continue;
576 std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>
582 double t = 0.00548952;
586 double p2 = energy * energy - m * m;
591 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
600 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
607 double bSelect = flat.fire();
608 if (bSelect < 0.819) {
610 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
613 else if (bSelect < 0.9954) {
615 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
622 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
629 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
634 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
646 std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom =
647 std::make_tuple(pdgid, m,
dirCalc(p, m));
648 v_prods.push_back(partMassMom);
652 for (
auto prodEntry : v_prods) {
655 ti_PDGID pdgid = std::get<0>(prodEntry);
656 td_Mass m = std::get<1>(prodEntry);
657 TLorentzVector pvec = std::get<2>(prodEntry);
659 std::string primary(
"primary");
663 if (pdgid == 1000020040) {
682 double costheta = (2.0 * flat.fire() - 1.0);
683 if (costheta < -1.0) costheta = -1.0;
684 if (costheta > 1.0) costheta = 1.0;
685 double const sintheta = sqrt(1.0 - costheta * costheta);
686 double const phi = 2.0 * M_PI * flat.fire();
687 return TLorentzVector{p * sintheta * std::cos(phi),
688 p * sintheta * std::sin(phi),
690 std::sqrt(p * p + m * m)};
696 std::string
const& filename,
697 std::string
const& PathToFile)
700 std::regex
const re_argon{
"42Ar.*"};
701 for (
size_t i = 0; i <
fNuclide.size(); i++) {
708 else if (std::regex_match(nuclide, re_argon) &&
fNuclide[i] ==
"42Ar") {
715 Bool_t addStatus = TH1::AddDirectoryStatus();
721 cet::search_path sp(
"FW_SEARCH_PATH");
722 std::string fn2 = PathToFile;
724 std::string fullname;
725 sp.find_file(fn2, fullname);
727 if (fullname.empty() || stat(fullname.c_str(), &sb) != 0)
729 <<
"Input spectrum file " << fn2 <<
" not found in FW_SEARCH_PATH!\n";
731 TFile
f(fullname.c_str(),
"READ");
732 TGraph* alphagraph = (TGraph*)
f.Get(
"Alphas");
733 TGraph* betagraph = (TGraph*)
f.Get(
"Betas");
734 TGraph* gammagraph = (TGraph*)
f.Get(
"Gammas");
735 TGraph* neutrongraph = (TGraph*)
f.Get(
"Neutrons");
738 int np = alphagraph->GetN();
739 double*
y = alphagraph->GetY();
744 auto alphahist = std::make_unique<TH1D>(name.c_str(),
"Alpha Spectrum", np, 0, np);
745 for (
int i = 0; i < np; i++) {
746 alphahist->SetBinContent(i + 1, y[i]);
747 alphahist->SetBinError(i + 1, 0);
758 int np = betagraph->GetN();
760 double*
y = betagraph->GetY();
765 auto betahist = std::make_unique<TH1D>(name.c_str(),
"Beta Spectrum", np, 0, np);
766 for (
int i = 0; i < np; i++) {
767 betahist->SetBinContent(i + 1, y[i]);
768 betahist->SetBinError(i + 1, 0);
779 int np = gammagraph->GetN();
780 double*
y = gammagraph->GetY();
785 auto gammahist = std::make_unique<TH1D>(name.c_str(),
"Gamma Spectrum", np, 0, np);
786 for (
int i = 0; i < np; i++) {
787 gammahist->SetBinContent(i + 1, y[i]);
788 gammahist->SetBinError(i + 1, 0);
799 int np = neutrongraph->GetN();
800 double*
y = neutrongraph->GetY();
805 auto neutronhist = std::make_unique<TH1D>(name.c_str(),
"Neutron Spectrum", np, 0, np);
806 for (
int i = 0; i < np; i++) {
807 neutronhist->SetBinContent(i + 1, y[i]);
808 neutronhist->SetBinError(i + 1, 0);
819 TH1::AddDirectory(addStatus);
845 throw cet::exception(
"RadioGen") <<
"Ununderstood nuclide: " << nuclide <<
"\n";
848 double rtype = flat.fire();
854 int itry = 0; itry < 10;
878 if (itype >= 0)
break;
882 <<
"Normalization problem with nuclide: " << nuclide <<
"\n";
886 if (p >= 0) { p = TMath::Sqrt(p); }
899 int nbinsx = hist.GetNbinsX();
900 std::vector<double> partialsum;
901 partialsum.resize(nbinsx + 1);
904 for (
int i = 1; i <=
nbinsx; i++) {
905 double hc = hist.GetBinContent(i);
907 throw cet::exception(
"RadioGen") <<
"Negative bin: " << i <<
" " << hist.GetName() <<
"\n";
908 partialsum[i] = partialsum[i - 1] + hc;
910 double integral = partialsum[
nbinsx];
911 if (integral == 0)
return 0;
913 for (
int i = 1; i <=
nbinsx; i++)
914 partialsum[i] /= integral;
916 double r1 = flat.fire();
917 int ibin = TMath::BinarySearch(nbinsx, &(partialsum[0]), r1);
918 Double_t
x = hist.GetBinLowEdge(ibin + 1);
919 if (r1 > partialsum[ibin]) {
920 x += hist.GetBinWidth(ibin + 1) * (r1 - partialsum[ibin]) /
921 (partialsum[ibin + 1] - partialsum[ibin]);
937 std::vector<double> vd_p = {.0015246};
938 for (
auto p : vd_p) {
939 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
947 std::vector<double> vd_p = {.0003126};
948 for (
auto p : vd_p) {
949 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
959 double chan1 = (0.052 / (0.052 + 0.020));
960 if (flat.fire() < chan1) {
961 std::vector<double> vd_p = {.0008997};
962 for (
auto p : vd_p) {
963 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
968 std::vector<double> vd_p = {.0024243};
969 for (
auto p : vd_p) {
970 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
980 double chan1 = (0.0033 / (0.0033 + 0.0201 + 0.041));
981 double chan2 = (0.0201 / (0.0033 + 0.0201 + 0.041));
982 double chanPick = flat.fire();
983 if (chanPick < chan1) {
984 std::vector<double> vd_p = {0.000692, 0.001228};
985 for (
auto p : vd_p) {
986 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
990 else if (chanPick < (chan1 + chan2)) {
991 std::vector<double> vd_p = {0.0010212};
992 for (
auto p : vd_p) {
993 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
998 std::vector<double> vd_p = {0.0019208};
999 for (
auto p : vd_p) {
1000 v_prods.emplace_back(pdgid, m,
dirCalc(p, m));
1028 template <
typename Stream>
1031 std::string
const& volumeName )
const 1034 out <<
"[#" << iNucl <<
"] " <<
fNuclide[iNucl] <<
" (" <<
fBq[iNucl] <<
" Bq/cm^3)" 1035 <<
" in " <<
fMaterial[iNucl] <<
" from " <<
fT0[iNucl] <<
" to " <<
fT1[iNucl]
1036 <<
" ns in ( " <<
fX0[iNucl] <<
", " <<
fY0[iNucl] <<
", " <<
fZ0[iNucl] <<
") to ( " 1037 <<
fX1[iNucl] <<
", " <<
fY1[iNucl] <<
", " <<
fZ1[iNucl] <<
")";
1038 if (!volumeName.empty()) out <<
" (\"" << volumeName <<
"\")";
1059 auto const trigTimeTick = timings.toTick<
electronics_tick>(timings.TriggerTime());
1063 return {double(timings.toTimeScale<
simulation_time>(trigTimeTick + beforeTicks)),
1064 double(timings.toTimeScale<
simulation_time>(trigTimeTick + afterTicks))};
code to link reconstructed objects back to the MC truth information
base_engine_t & createEngine(seed_t seed)
std::vector< std::unique_ptr< TH1D > > gammaspectrum
TLorentzVector dirCalc(double p, double m)
Utilities related to art service access.
timescale_traits< ElectronicsTimeCategory >::tick_t electronics_tick
A point on the electronics time scale expressed in its ticks.
bool apply(geo::GeoNodePath &path, Op &&op) const
Applies the specified operation to all nodes under the path.
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
std::vector< double > fZ0
Bottom corner z position (cm) in world coordinates.
void SetOrigin(simb::Origin_t origin)
Executes an operation on all the nodes of the ROOT geometry.
Definition of util::enumerate().
std::vector< std::string > fNuclide
List of nuclides to simulate. Example: "39Ar".
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::pair< double, double > defaulttimewindow() const
Returns the start and end of the readout window.
Class representing a path in ROOT geometry.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
unsigned int ReadOutWindowSize() const
Class representing a path in ROOT geometry.
void Ar42Gamma5(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
void dumpNuclideSettings(Stream &&out, std::size_t iNucl, std::string const &volumeName={}) const
Prints the settings for the specified nuclide and volume.
TGeoManager * ROOTGeoManager() const
Access to the ROOT geometry description manager.
std::vector< double > neutronintegral
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
std::vector< std::string > spectrumname
Interface to detinfo::DetectorClocks.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Access the description of detector geometry.
std::vector< std::unique_ptr< TH1D > > alphaspectrum
std::vector< double > fT1
End of time window to simulate in ns.
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
std::vector< double > fY1
Top corner y position (cm) in world coordinates.
timescale_traits< ElectronicsTimeCategory >::tick_interval_t electronics_time_ticks
An interval on the electronics time scale expressed in its ticks.
#define DEFINE_ART_MODULE(klass)
Definitions of geometry vector data types.
timescale_traits< SimulationTimeCategory >::time_point_t simulation_time
A point in time on the simulation time scale.
std::vector< double > fZ1
Top corner z position (cm) in world coordinates.
single particles thrown at the detector
detinfo::DetectorTimings makeDetectorTimings(detinfo::DetectorClocksData const &detClocks)
int trackidcounter
Serial number for the MC track ID.
Utilities to extend the interface of geometry vectors.
void SampleOne(unsigned int i, simb::MCTruth &mct)
Checks that the node represents a box well aligned to world frame axes.
void Ar42Gamma3(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
unsigned int NumberTimeSamples() const
void readfile(std::string nuclide, std::string const &filename, std::string const &PathToFile)
Test of util::counter and support utilities.
void beginRun(art::Run &run)
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
Module to generate particles created by radiological decay, patterned off of SingleGen.
double samplefromth1d(TH1D &hist)
std::vector< double > alphaintegral
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
std::size_t addvolume(std::string const &volumeName)
void samplespectrum(std::string nuclide, int &itype, double &t, double &m, double &p)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
std::vector< double > betaintegral
void Add(simb::MCParticle const &part)
void Ar42Gamma2(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
CLHEP::HepRandomEngine & fEngine
std::vector< double > fT0
Beginning of time window to simulate in ns.
Representation of a node and its ancestry.
void Ar42Gamma4(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
std::vector< std::unique_ptr< TH1D > > neutronspectrum
std::vector< std::string > fMaterial
List of regexes of materials in which to generate the decays. Example: "LAr".
std::string const & DetectorName() const
Returns a string with the name of the detector, as configured.
std::vector< double > fX0
Bottom corner x position (cm) in world coordinates.
Event generator information.
bool fIsFirstSignalSpecial
std::vector< double > gammaintegral
std::vector< double > fBq
Radioactivity in Becquerels (decay per sec) per cubic cm.
Namespace collecting geometry-related classes utilities.
Event Generation using GENIE, cosmics or single particles.
Namespace including different time scales as defined in LArSoft.
std::vector< std::unique_ptr< TH1D > > betaspectrum
std::vector< double > fX1
Top corner x position (cm) in world coordinates.
art framework interface to geometry description
ROOT::Math::Transform3D TransformationMatrix
Type of transformation matrix used in geometry.
constexpr Point origin()
Returns a origin position with a point of the specified type.
cet::coded_exception< error, detail::translate > exception
std::vector< double > fY0
Bottom corner y position (cm) in world coordinates.
void produce(art::Event &evt)