LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
RadioGen_module.cc
Go to the documentation of this file.
1 // Rn222 generation feature added by gleb.sinev@duke.edu
9 // (based on a generator by jason.stock@mines.sdsmt.edu)
10 // JStock. Added preliminary changes to get ready for Ar42 implimentation. This includes allowing for multiple particles from the same decay.
11 // Ar42 generation added by JStock (jason.stock@mines.sdsmt.edu).
12 // Ar42 is designed to handle 5 separate different
13 // beta decay modes, each with it's own chains of
14 // possible dexcitation gammas. Because of the
15 // high energies, and the relatively high rate
16 // expected in the DUNE FD, these chains are
17 // completely simulated instead of relying on the
18 // existing machinery in this module. To make the
19 // treatment of multiple decay products and
20 // dexcitation chains generally available would
21 // require a significant redesign of the module
22 // and possibly of the .root files data structure
23 // for each radiological.
24 //
25 // Dec 01, 2017 JStock
26 // Adding the ability to make 8.997 MeV Gammas for
27 // Ni59 Calibration sources. This is another "hacky"
28 // fix to something that really deserves a more
29 // elegant and comprehensive solution.
30 //
31 // Aug 2023 L Paulucci
32 // Adding a fhicl parameter for generic path to radiogen files
34 
35 // C++ includes.
36 #include <cassert>
37 #include <cmath>
38 #include <iterator>
39 #include <memory>
40 #include <regex>
41 #include <string>
42 #include <sys/stat.h>
43 #include <utility> // std::pair<>, std::tuple<>
44 #include <vector>
45 
46 // Framework includes
53 #include "cetlib/exempt_ptr.h"
54 #include "cetlib/search_path.h"
55 #include "cetlib_except/exception.h"
56 #include "fhiclcpp/ParameterSet.h"
58 
59 // nurandom includes
61 
62 // nusimdata, nugen includes
66 
67 // lar includes
76 #include "larcorealg/Geometry/geo_vectors_utils.h" // geo::...::makeFromCoords()
85 
86 // root includes
87 
88 #include "TFile.h"
89 #include "TGenPhaseSpace.h"
90 #include "TGeoBBox.h"
91 #include "TGeoManager.h"
92 #include "TGeoMaterial.h"
93 #include "TGeoNode.h"
94 #include "TGeoVolume.h"
95 #include "TGraph.h"
96 #include "TH1.h"
97 #include "TH1D.h"
98 #include "TLorentzVector.h"
99 #include "TMath.h"
100 
101 #include "CLHEP/Random/RandFlat.h"
102 #include "CLHEP/Random/RandPoisson.h"
103 
104 namespace simb {
105  class MCTruth;
106 }
107 
108 namespace evgen {
109 
193  class RadioGen : public art::EDProducer {
194  public:
195  explicit RadioGen(fhicl::ParameterSet const& pset);
196 
197  private:
198  // This is called for each event.
199  void produce(art::Event& evt);
200  void beginRun(art::Run& run);
201 
202  typedef int
203  ti_PDGID; // These typedefs may look odd, and unecessary. I chose to use them to make the tuples I use later more readable. ti, type integer :JStock
204  typedef double
205  td_Mass; // These typedefs may look odd, and unecessary. I chose to use them to make the tuples I use later more readable. td, type double :JStock
206 
209  std::size_t addvolume(std::string const& volumeName);
210 
212  void SampleOne(unsigned int i, simb::MCTruth& mct);
213 
214  TLorentzVector dirCalc(double p, double m);
215 
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);
218 
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);
223 
224  // recoded so as to use the LArSoft-managed random number generator
225  double samplefromth1d(TH1D& hist);
226 
228  template <typename Stream>
229  void dumpNuclideSettings(Stream&& out,
230  std::size_t iNucl,
231  std::string const& volumeName = {}) const;
232 
234  std::pair<double, double> defaulttimewindow() const;
235 
236  // itype = pdg code: 1000020040: alpha. itype=11: beta. -11: positron, itype=22: gamma. -1: error
237  // t is the kinetic energy in GeV (= E-m)
238  // m = mass in GeV
239  // p = momentum in GeV
240 
241  //double betaphasespace(double mass, double q); // older parameterization.
242 
243  // the generator randomly samples points in a rectangular prism of space and time, and only selects those points in
244  // volumes with materials that match the regexes in fMaterial. One can use wildcards * and ? for broader matches.
245 
246  std::vector<std::string> fNuclide;
247  std::vector<std::string>
249  std::vector<double> fBq;
250  std::vector<double> fT0;
251  std::vector<double> fT1;
252  std::vector<double> fX0;
253  std::vector<double> fY0;
254  std::vector<double> fZ0;
255  std::vector<double> fX1;
256  std::vector<double> fY1;
257  std::vector<double> fZ1;
260  std::string fPathToFile;
261 
262  // leftovers from the phase space generator
263  // const double gevperamu = 0.931494061;
264  // TGenPhaseSpace rg; // put this here so we don't constantly construct and destruct it
265 
266  std::vector<std::string> spectrumname;
267  std::vector<std::unique_ptr<TH1D>> alphaspectrum;
268  std::vector<double> alphaintegral;
269  std::vector<std::unique_ptr<TH1D>> betaspectrum;
270  std::vector<double> betaintegral;
271  std::vector<std::unique_ptr<TH1D>> gammaspectrum;
272  std::vector<double> gammaintegral;
273  std::vector<std::unique_ptr<TH1D>> neutronspectrum;
274  std::vector<double> neutronintegral;
275  CLHEP::HepRandomEngine& fEngine;
276  };
277 }
278 
279 namespace {
280  constexpr double m_e = 0.000510998928; // mass of electron in GeV
281  constexpr double m_alpha = 3.727379240; // mass of an alpha particle in GeV
282  constexpr double m_neutron = 0.9395654133; // mass of a neutron in GeV
283 }
284 
285 namespace evgen {
286 
287  //____________________________________________________________________________
288  RadioGen::RadioGen(fhicl::ParameterSet const& pset)
289  : EDProducer{pset}
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", {})}
296  , fIsFirstSignalSpecial{pset.get<bool>("IsFirstSignalSpecial", false)}
297  // create a default random engine; obtain the random seed from NuRandomService,
298  // unless overridden in configuration with key "Seed"
299  , fEngine(art::ServiceHandle<rndm::NuRandomService>()->registerAndSeedEngine(createEngine(0),
300  pset,
301  "Seed"))
302  {
303  produces<std::vector<simb::MCTruth>>();
304  produces<sumdata::RunData, art::InRun>();
305 
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", {});
312 
313  if (fPathToFile.find_last_of('/') !=
314  fPathToFile.length() - 1) { //Checking that last character in path is a forward slash
315  throw cet::exception("RadioGen")
316  << "Last entry of the path to radiological files needs to be a forward slash ('/')\n";
317  }
318 
319  if (fT0.empty() || fT1.empty()) { // better be both empty...
320  if (!fT0.empty() || !fT1.empty()) {
322  << "RadioGen T0 and T1 need to be both non-empty, or both empty"
323  " (now T0 has "
324  << fT0.size() << " entries and T1 has " << fT0.size() << ")\n";
325  }
326  auto const [defaultT0, defaultT1] = defaulttimewindow();
327  t0.push_back(defaultT0);
328  t1.push_back(defaultT1);
329  }
330 
331  //
332  // First, the materials are assigned to the coordinates explicitly
333  // configured; then, each of the remaining materials is assigned to one
334  // of the named volumes.
335  // TODO: reaction to mismatches is "not so great".
336  //
337  mf::LogInfo("RadioGen") << "Configuring activity:";
338  for (std::size_t iCoord : util::counter(fX0.size())) {
339 
340  fNuclide.push_back(nuclide.at(iCoord));
341  fMaterial.push_back(material.at(iCoord));
342  fBq.push_back(Bq.at(iCoord));
343  // replicate the last timing if none is specified
344  fT0.push_back(t0[std::min(iCoord, t0.size() - 1U)]);
345  fT1.push_back(t1[std::min(iCoord, t1.size() - 1U)]);
346 
347  dumpNuclideSettings(mf::LogVerbatim("RadioGen"), iCoord);
348 
349  } // for
350 
351  std::size_t iNext = fX0.size();
352  auto const volumes = pset.get<std::vector<std::string>>("Volumes", {});
353  for (auto&& [iVolume, volName] : util::enumerate(volumes)) {
354  // this expands the coordinate vectors
355  auto const nVolumes = addvolume(volName);
356  if (nVolumes == 0) {
358  << "No volume named '" << volName << "' was found!\n";
359  }
360 
361  std::size_t const iVolFirst = fNuclide.size();
362  for (auto iVolInstance : util::counter(nVolumes)) {
363  fNuclide.push_back(nuclide.at(iNext));
364  fMaterial.push_back(material.at(iNext));
365  fBq.push_back(Bq.at(iNext));
366  // replicate the last timing if none is specified
367  fT0.push_back(t0[std::min(iNext, t0.size() - 1U)]);
368  fT1.push_back(t1[std::min(iNext, t1.size() - 1U)]);
369 
370  dumpNuclideSettings(mf::LogVerbatim("RadioGen"), iVolFirst + iVolInstance, volName);
371  }
372  ++iNext;
373  } // for
374 
375  // check for consistency of vector sizes
376  unsigned int nsize = fNuclide.size();
377  if (nsize == 0) {
378  throw art::Exception(art::errors::Configuration) << "No nuclide configured!\n";
379  }
380 
381  if (fMaterial.size() != nsize)
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";
401 
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);
406  }
407  else if (nuclideName == "85Kr") {
408  readfile("85Kr", "Krypton_85.root", fPathToFile);
409  }
410  else if (nuclideName == "40K") {
411  readfile("40K", "Potassium_40.root", fPathToFile);
412  }
413  else if (nuclideName == "232Th") {
414  readfile("232Th", "Thorium_232.root", fPathToFile);
415  }
416  else if (nuclideName == "238U") {
417  readfile("238U", "Uranium_238.root", fPathToFile);
418  }
419  else if (nuclideName == "222Rn") {
420  continue;
421  } //Rn222 is handled separately later
422  else if (nuclideName == "59Ni") {
423  continue;
424  } //Rn222 is handled separately later
425  else if (nuclideName == "42Ar") {
426  readfile(
427  "42Ar_1",
428  "Argon_42_1.root",
429  fPathToFile); //Each possible beta decay mode of Ar42 is given it's own .root file for now.
430  readfile(
431  "42Ar_2",
432  "Argon_42_2.root",
433  fPathToFile); //This allows us to know which decay chain to follow for the dexcitation gammas.
434  readfile(
435  "42Ar_3",
436  "Argon_42_3.root",
437  fPathToFile); //The dexcitation gammas are not included in the root files as we want to
438  readfile(
439  "42Ar_4",
440  "Argon_42_4.root",
441  fPathToFile); //probabilistically simulate the correct coincident gammas, which we cannot guarantee
442  readfile("42Ar_5", "Argon_42_5.root", fPathToFile); //by sampling a histogram.
443  continue;
444  } //Ar42 is handeled separately later
445  else {
446  std::string searchName = nuclideName;
447  searchName += ".root";
448  readfile(nuclideName, searchName, fPathToFile);
449  }
450  }
451  }
452 
453  //____________________________________________________________________________
455  {
457  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
458  }
459 
460  //____________________________________________________________________________
462  {
464  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
465 
466  simb::MCTruth truth;
468 
469  trackidcounter = -1;
470  for (unsigned int i = 0; i < fNuclide.size(); ++i) {
471  SampleOne(i, truth);
472  } //end loop over nuclides
473 
474  MF_LOG_DEBUG("RadioGen") << truth;
475  truthcol->push_back(truth);
476  evt.put(std::move(truthcol));
477  }
478 
479  //____________________________________________________________________________
480  std::size_t RadioGen::addvolume(std::string const& volumeName)
481  {
482 
483  auto const& geom = *(lar::providerFrom<geo::Geometry>());
484 
485  std::vector<geo::GeoNodePath> volumePaths;
486  auto findVolume = [&volumePaths, volumeName](auto& path) {
487  if (path.current().GetVolume()->GetName() == volumeName) volumePaths.push_back(path);
488  return true;
489  };
490 
491  geo::ROOTGeometryNavigator navigator{*(geom.ROOTGeoManager())};
492 
493  navigator.apply(findVolume);
494 
495  for (geo::GeoNodePath const& path : volumePaths) {
496 
497  //
498  // find the coordinates of the volume in local coordinates
499  //
500  TGeoShape const* pShape = path.current().GetVolume()->GetShape();
501  auto pBox = dynamic_cast<TGeoBBox const*>(pShape);
502  if (!pBox) {
503  throw cet::exception("RadioGen") << "Volume '" << path.current().GetName() << "' is a "
504  << pShape->IsA()->GetName() << ", not a TGeoBBox.\n";
505  }
506 
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()};
509 
510  //
511  // convert to world coordinates
512  //
513 
514  auto const trans = path.currentTransformation<geo::TransformationMatrix>();
515 
516  geo::Point_t min, max;
517  trans.Transform(origin - diag, min);
518  trans.Transform(origin + diag, max);
519 
520  //
521  // add to the coordinates
522  //
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()));
529 
530  } // for
531 
532  return volumePaths.size();
533  } // RadioGen::addvolume()
534 
535  //____________________________________________________________________________
536  // Generate radioactive decays per nuclide per volume according to the FCL parameters
537 
538  void RadioGen::SampleOne(unsigned int i, simb::MCTruth& mct)
539  {
541  TGeoManager* geomanager = geo->ROOTGeoManager();
542 
543  CLHEP::RandFlat flat(fEngine);
544  CLHEP::RandPoisson poisson(fEngine);
545 
546  // figure out how many decays to generate, assuming that the entire prism consists of the radioactive material.
547  // we will skip over decays in other materials later.
548 
549  double rate =
550  fabs(fBq[i] * (fT1[i] - fT0[i]) * (fX1[i] - fX0[i]) * (fY1[i] - fY0[i]) * (fZ1[i] - fZ0[i])) /
551  1.0E9;
552  long ndecays = poisson.shoot(rate);
553 
554  std::regex const re_material{fMaterial[i]};
555  for (unsigned int idecay = 0; idecay < ndecays; idecay++) {
556  // generate just one particle at a time. Need a little recoding if a radioactive
557  // decay generates multiple daughters that need simulation
558  // uniformly distributed in position and time
559  //
560  // JStock: Leaving this as a single position for the decay products. For now I will assume they all come from the same spot.
561  TLorentzVector pos(
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]),
565  (idecay == 0 && fIsFirstSignalSpecial) ? 0 : (fT0[i] + flat.fire() * (fT1[i] - fT0[i])));
566 
567  // discard decays that are not in the proper material
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;
571 
572  //Moved pdgid into the next statement, so that it is localized.
573  // electron=11, photon=22, alpha = 1000020040, neutron = 2112
574 
575  //JStock: Allow us to have different particles from the same decay. This requires multiple momenta.
576  std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>
577  v_prods; //(First is for PDGID, second is mass, third is Momentum)
578 
579  if (fNuclide[i] == "222Rn") // Treat 222Rn separately
580  {
581  double p = 0;
582  double t = 0.00548952;
583  td_Mass m = m_alpha;
584  ti_PDGID pdgid = 1000020040; //td_Mass = double. ti_PDGID = int;
585  double energy = t + m;
586  double p2 = energy * energy - m * m;
587  if (p2 > 0)
588  p = TMath::Sqrt(p2);
589  else
590  p = 0;
591  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
592  } //End special case RN222
593  else if (
594  fNuclide[i] ==
595  "59Ni") { //Treat 59Ni Calibration Source separately (as I haven't made a spectrum for it, and ultimately it should be handeled with multiple particle outputs.
596  double p = 0.008997;
597  td_Mass m = 0;
598  ti_PDGID pdgid =
599  22; // td_Mas=double. ti_PDFID=int. Assigning p directly, as t=p for gammas.
600  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
601  } //end special case Ni59 calibration source
602  else if (fNuclide[i] == "42Ar") { // Spot for special treatment of Ar42.
603  double p = 0;
604  double t = 0;
605  td_Mass m = 0;
606  ti_PDGID pdgid = 0; //td_Mass = double. ti_PDGID = int;
607  double bSelect = flat.fire(); //Make this a random number from 0 to 1.
608  if (bSelect < 0.819) { //beta channel 1. No Gamma. beta Q value 3525.22 keV
609  samplespectrum("42Ar_1", pdgid, t, m, p);
610  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
611  //No gamma here.
612  }
613  else if (bSelect < 0.9954) { //beta channel 2. 1 Gamma (1524.6 keV). beta Q value 2000.62
614  samplespectrum("42Ar_2", pdgid, t, m, p);
615  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
616  Ar42Gamma2(v_prods);
617  }
618  else if (
619  bSelect <
620  0.9988) { //beta channel 3. 1 Gamma Channel. 312.6 keV + gamma 2. beta Q value 1688.02 keV
621  samplespectrum("42Ar_3", pdgid, t, m, p);
622  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
623  Ar42Gamma3(v_prods);
624  }
625  else if (
626  bSelect <
627  0.9993) { //beta channel 4. 2 Gamma Channels. Either 899.7 keV (i 0.052) + gamma 2 or 2424.3 keV (i 0.020). beta Q value 1100.92 keV
628  samplespectrum("42Ar_4", pdgid, t, m, p);
629  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
630  Ar42Gamma4(v_prods);
631  }
632  else { //beta channel 5. 3 gamma channels. 692.0 keV + 1228.0 keV + Gamma 2 (i 0.0033) ||OR|| 1021.2 keV + gamma 4 (i 0.0201) ||OR|| 1920.8 keV + gamma 2 (i 0.041). beta Q value 79.82 keV
633  samplespectrum("42Ar_5", pdgid, t, m, p);
634  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
635  Ar42Gamma5(v_prods);
636  }
637  //Add beta.
638  //Call gamma function for beta mode.
639  }
640  else { //General Case.
641  double p = 0;
642  double t = 0;
643  td_Mass m = 0;
644  ti_PDGID pdgid = 0; //td_Mass = double. ti_PDGID = int;
645  samplespectrum(fNuclide[i], pdgid, t, m, p);
646  std::tuple<ti_PDGID, td_Mass, TLorentzVector> partMassMom =
647  std::make_tuple(pdgid, m, dirCalc(p, m));
648  v_prods.push_back(partMassMom);
649  } //end else (not RN or other special case
650 
651  //JStock: Modify this to now loop over the v_prods.
652  for (auto prodEntry : v_prods) {
653  // set track id to a negative serial number as these are all primary particles and have id <= 0
654  int trackid = trackidcounter;
655  ti_PDGID pdgid = std::get<0>(prodEntry);
656  td_Mass m = std::get<1>(prodEntry);
657  TLorentzVector pvec = std::get<2>(prodEntry);
658  trackidcounter--;
659  std::string primary("primary");
660 
661  // alpha particles need a little help since they're not in the TDatabasePDG table
662  // // so don't rely so heavily on default arguments to the MCParticle constructor
663  if (pdgid == 1000020040) {
664  simb::MCParticle part(trackid, pdgid, primary, -1, m, 1);
665  part.AddTrajectoryPoint(pos, pvec);
666  mct.Add(part);
667  } // end "If alpha"
668  else {
669  simb::MCParticle part(trackid, pdgid, primary);
670  part.AddTrajectoryPoint(pos, pvec);
671  mct.Add(part);
672  } // end All standard cases.
673  } //End Loop over all particles produces in this single decay.
674  }
675  }
676 
677  //Calculate an arbitrary direction with a given magnitude p
678  TLorentzVector RadioGen::dirCalc(double p, double m)
679  {
680  CLHEP::RandFlat flat(fEngine);
681  // isotropic production angle for the decay product
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),
689  p * costheta,
690  std::sqrt(p * p + m * m)};
691  }
692 
693  // only reads those files that are on the fNuclide list. Copy information from the TGraphs to TH1D's
694 
695  void RadioGen::readfile(std::string nuclide,
696  std::string const& filename,
697  std::string const& PathToFile)
698  {
699  bool found{false};
700  std::regex const re_argon{"42Ar.*"};
701  for (size_t i = 0; i < fNuclide.size(); i++) {
702  if (
703  fNuclide[i] ==
704  nuclide) { //This check makes sure that the nuclide we are searching for is in fact in our fNuclide list. Ar42 handeled separately.
705  found = true;
706  break;
707  } //End If nuclide is in our list. Next is the special case of Ar42
708  else if (std::regex_match(nuclide, re_argon) && fNuclide[i] == "42Ar") {
709  found = true;
710  break;
711  }
712  }
713  if (!found) return;
714 
715  Bool_t addStatus = TH1::AddDirectoryStatus();
716  TH1::AddDirectory(
717  kFALSE); // cloned histograms go in memory, and aren't deleted when files are closed.
718  // be sure to restore this state when we're out of the routine.
719 
720  spectrumname.push_back(nuclide);
721  cet::search_path sp("FW_SEARCH_PATH");
722  std::string fn2 = PathToFile;
723  fn2 += filename;
724  std::string fullname;
725  sp.find_file(fn2, fullname);
726  struct stat sb;
727  if (fullname.empty() || stat(fullname.c_str(), &sb) != 0)
728  throw cet::exception("RadioGen")
729  << "Input spectrum file " << fn2 << " not found in FW_SEARCH_PATH!\n";
730 
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");
736 
737  if (alphagraph) {
738  int np = alphagraph->GetN();
739  double* y = alphagraph->GetY();
740  std::string name;
741  name = "RadioGen_";
742  name += nuclide;
743  name += "_Alpha";
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);
748  }
749  alphaintegral.push_back(alphahist->Integral());
750  alphaspectrum.push_back(move(alphahist));
751  }
752  else {
753  alphaintegral.push_back(0);
754  alphaspectrum.push_back(nullptr);
755  }
756 
757  if (betagraph) {
758  int np = betagraph->GetN();
759 
760  double* y = betagraph->GetY();
761  std::string name;
762  name = "RadioGen_";
763  name += nuclide;
764  name += "_Beta";
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);
769  }
770  betaintegral.push_back(betahist->Integral());
771  betaspectrum.push_back(move(betahist));
772  }
773  else {
774  betaintegral.push_back(0);
775  betaspectrum.push_back(nullptr);
776  }
777 
778  if (gammagraph) {
779  int np = gammagraph->GetN();
780  double* y = gammagraph->GetY();
781  std::string name;
782  name = "RadioGen_";
783  name += nuclide;
784  name += "_Gamma";
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);
789  }
790  gammaintegral.push_back(gammahist->Integral());
791  gammaspectrum.push_back(move(gammahist));
792  }
793  else {
794  gammaintegral.push_back(0);
795  gammaspectrum.push_back(nullptr);
796  }
797 
798  if (neutrongraph) {
799  int np = neutrongraph->GetN();
800  double* y = neutrongraph->GetY();
801  std::string name;
802  name = "RadioGen_";
803  name += nuclide;
804  name += "_Neutron";
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);
809  }
810  neutronintegral.push_back(neutronhist->Integral());
811  neutronspectrum.push_back(move(neutronhist));
812  }
813  else {
814  neutronintegral.push_back(0);
815  neutronspectrum.push_back(nullptr);
816  }
817 
818  f.Close();
819  TH1::AddDirectory(addStatus);
820 
821  double total =
822  alphaintegral.back() + betaintegral.back() + gammaintegral.back() + neutronintegral.back();
823  if (total > 0) {
824  alphaintegral.back() /= total;
825  betaintegral.back() /= total;
826  gammaintegral.back() /= total;
827  neutronintegral.back() /= total;
828  }
829  }
830 
831  void RadioGen::samplespectrum(std::string nuclide, int& itype, double& t, double& m, double& p)
832  {
833  CLHEP::RandFlat flat(fEngine);
834 
835  int inuc = -1;
836  for (size_t i = 0; i < spectrumname.size(); i++) {
837  if (nuclide == spectrumname[i]) {
838  inuc = i;
839  break;
840  }
841  }
842  if (inuc == -1) {
843  t = 0; // throw an exception in the future
844  itype = 0;
845  throw cet::exception("RadioGen") << "Ununderstood nuclide: " << nuclide << "\n";
846  }
847 
848  double rtype = flat.fire();
849 
850  itype = -1;
851  m = 0;
852  p = 0;
853  for (
854  int itry = 0; itry < 10;
855  itry++) // maybe a tiny normalization issue with a sum of 0.99999999999 or something, so try a few times.
856  {
857  if (rtype <= alphaintegral[inuc] && alphaspectrum[inuc] != nullptr) {
858  itype = 1000020040; // alpha
859  m = m_alpha;
860  t = samplefromth1d(*alphaspectrum[inuc]) / 1000000.0;
861  }
862  else if (rtype <= alphaintegral[inuc] + betaintegral[inuc] && betaspectrum[inuc] != nullptr) {
863  itype = 11; // beta
864  m = m_e;
865  t = samplefromth1d(*betaspectrum[inuc]) / 1000000.0;
866  }
867  else if (rtype <= alphaintegral[inuc] + betaintegral[inuc] + gammaintegral[inuc] &&
868  gammaspectrum[inuc] != nullptr) {
869  itype = 22; // gamma
870  m = 0;
871  t = samplefromth1d(*gammaspectrum[inuc]) / 1000000.0;
872  }
873  else if (neutronspectrum[inuc] != nullptr) {
874  itype = 2112;
875  m = m_neutron;
876  t = samplefromth1d(*neutronspectrum[inuc]) / 1000000.0;
877  }
878  if (itype >= 0) break;
879  }
880  if (itype == -1) {
881  throw cet::exception("RadioGen")
882  << "Normalization problem with nuclide: " << nuclide << "\n";
883  }
884  double e = t + m;
885  p = e * e - m * m;
886  if (p >= 0) { p = TMath::Sqrt(p); }
887  else {
888  p = 0;
889  }
890  }
891 
892  // this is just a copy of TH1::GetRandom that uses the art-managed CLHEP random number generator instead of gRandom
893  // and a better handling of negative bin contents
894 
896  {
897  CLHEP::RandFlat flat(fEngine);
898 
899  int nbinsx = hist.GetNbinsX();
900  std::vector<double> partialsum;
901  partialsum.resize(nbinsx + 1);
902  partialsum[0] = 0;
903 
904  for (int i = 1; i <= nbinsx; i++) {
905  double hc = hist.GetBinContent(i);
906  if (hc < 0)
907  throw cet::exception("RadioGen") << "Negative bin: " << i << " " << hist.GetName() << "\n";
908  partialsum[i] = partialsum[i - 1] + hc;
909  }
910  double integral = partialsum[nbinsx];
911  if (integral == 0) return 0;
912  // normalize to unit sum
913  for (int i = 1; i <= nbinsx; i++)
914  partialsum[i] /= integral;
915 
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]);
922  }
923  return x;
924  }
925 
926  //Ar42 uses BNL tables for K-42 from Aug 2017
927  //beta channel 1. No Gamma. beta Q value 3525.22 keV
928  //beta channel 2. 1 Gamma (1524.6 keV). beta Q value 2000.62
929  //beta channel 3. 1 Gamma Channel. 312.6 keV + gamma 2. beta Q value 1688.02 keV
930  //beta channel 4. 2 Gamma Channels. Either 899.7 keV (i 0.052) + gamma 2 or 2424.3 keV (i 0.020). beta Q value 1100.92 keV
931  //beta channel 5. 3 gamma channels. 692.0 keV + 1228.0 keV + Gamma 2 (i 0.0033) ||OR|| 1021.2 keV + gamma 4 (i 0.0201) ||OR|| 1920.8 keV + gamma 2 (i 0.041). beta Q value 79.82 keV
932  //No Ar42Gamma1 as beta channel 1 does not produce a dexcitation gamma.
933  void RadioGen::Ar42Gamma2(std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>& v_prods)
934  {
935  ti_PDGID pdgid = 22;
936  td_Mass m = 0.0; //we are writing gammas
937  std::vector<double> vd_p = {.0015246}; //Momentum in GeV
938  for (auto p : vd_p) {
939  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
940  }
941  }
942 
943  void RadioGen::Ar42Gamma3(std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>& v_prods)
944  {
945  ti_PDGID pdgid = 22;
946  td_Mass m = 0.0; //we are writing gammas
947  std::vector<double> vd_p = {.0003126};
948  for (auto p : vd_p) {
949  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
950  }
951  Ar42Gamma2(v_prods);
952  }
953 
954  void RadioGen::Ar42Gamma4(std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>& v_prods)
955  {
956  CLHEP::RandFlat flat(fEngine);
957  ti_PDGID pdgid = 22;
958  td_Mass m = 0.0; //we are writing gammas
959  double chan1 = (0.052 / (0.052 + 0.020));
960  if (flat.fire() < chan1) {
961  std::vector<double> vd_p = {.0008997}; //Momentum in GeV
962  for (auto p : vd_p) {
963  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
964  }
965  Ar42Gamma2(v_prods);
966  }
967  else {
968  std::vector<double> vd_p = {.0024243}; //Momentum in GeV
969  for (auto p : vd_p) {
970  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
971  }
972  }
973  }
974 
975  void RadioGen::Ar42Gamma5(std::vector<std::tuple<ti_PDGID, td_Mass, TLorentzVector>>& v_prods)
976  {
977  CLHEP::RandFlat flat(fEngine);
978  ti_PDGID pdgid = 22;
979  td_Mass m = 0.0; //we are writing gammas
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}; //Momentum in GeV
985  for (auto p : vd_p) {
986  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
987  }
988  Ar42Gamma2(v_prods);
989  }
990  else if (chanPick < (chan1 + chan2)) {
991  std::vector<double> vd_p = {0.0010212}; //Momentum in GeV
992  for (auto p : vd_p) {
993  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
994  }
995  Ar42Gamma4(v_prods);
996  }
997  else {
998  std::vector<double> vd_p = {0.0019208}; //Momentum in GeV
999  for (auto p : vd_p) {
1000  v_prods.emplace_back(pdgid, m, dirCalc(p, m));
1001  }
1002  Ar42Gamma2(v_prods);
1003  }
1004  }
1005 
1006  // phase space generator for beta decay -- keep it as a comment in case we ever want to revive it
1007 
1008  // double RadioGen::betaphasespace(double mass, double q)
1009  //{
1010  // CLHEP::RandFlat flat(fEngine);
1011  // double p = 0;
1012  // double mi = mass+q+m_e;
1013  // TLorentzVector p0(0,0,0,mi); // pre-decay four-vector
1014  // double masses[3] = {0,m_e,mass}; // neutrino, electron, nucleus
1015  // rg.SetDecay(p0,3,masses);
1016  // double wmax = rg.GetWtMax();
1017  // for (int igen=0;igen<1000;igen++) // cap the retries at 1000
1018  // {
1019  // double weight = rg.Generate(); // want to unweight these if possible
1020  // TLorentzVector *e4v = rg.GetDecay(1); // get electron four-vector
1021  // p = e4v->P();
1022  // if (weight >= wmax * flat.fire()) break;
1023  // }
1024  //return p;
1025  //}
1026 
1027  //____________________________________________________________________________
1028  template <typename Stream>
1030  std::size_t iNucl,
1031  std::string const& volumeName /* = {} */) const
1032  {
1033 
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 << "\")";
1039 
1040  } // RadioGen::dumpNuclideSettings()
1041 
1042  //____________________________________________________________________________
1043  std::pair<double, double> RadioGen::defaulttimewindow() const
1044  {
1045  // values are in simulation time scale (and therefore in [ns])
1046  using namespace detinfo::timescales;
1047 
1048  auto const& timings = detinfo::makeDetectorTimings(
1050  detinfo::DetectorPropertiesData const& detInfo =
1052 
1053  //
1054  // we take a number of (TPC electronics) ticks before the trigger time,
1055  // and we go to a number of ticks after the trigger time;
1056  // that shift is one readout window size
1057  //
1058 
1059  auto const trigTimeTick = timings.toTick<electronics_tick>(timings.TriggerTime());
1060  electronics_time_ticks const beforeTicks{-static_cast<int>(detInfo.ReadOutWindowSize())};
1061  electronics_time_ticks const afterTicks{detInfo.NumberTimeSamples()};
1062 
1063  return {double(timings.toTimeScale<simulation_time>(trigTimeTick + beforeTicks)),
1064  double(timings.toTimeScale<simulation_time>(trigTimeTick + afterTicks))};
1065  } // RadioGen::defaulttimewindow()
1066 
1067 } //end namespace evgen
1068 
Float_t x
Definition: compare.C:6
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.
TTree * t1
Definition: plottest35.C:26
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)
Definition: MCParticle.h:258
constexpr auto fullRun()
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
std::vector< double > fZ0
Bottom corner z position (cm) in world coordinates.
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
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.
Float_t y
Definition: compare.C:6
Class representing a path in ROOT geometry.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Run.h:121
Class representing a path in ROOT geometry.
void Ar42Gamma5(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
Particle class.
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
Definition: Run.h:37
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
Definition: enumerate.h:65
std::vector< std::string > spectrumname
Interface to detinfo::DetectorClocks.
TFile f
Definition: plotHisto.C:6
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
Access the description of detector geometry.
TString part[npart]
Definition: Style.C:32
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.
Definition: counter.h:295
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)
Definition: ModuleMacros.h:65
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
Definition: MCTruth.h:26
detinfo::DetectorTimings makeDetectorTimings(detinfo::DetectorClocksData const &detClocks)
int trackidcounter
Serial number for the MC track ID.
double energy
Definition: plottest35.C:25
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)
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.
ART objects.
double samplefromth1d(TH1D &hist)
std::vector< double > alphaintegral
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
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.
Definition: geo_vectors.h:180
std::vector< double > betaintegral
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
void Ar42Gamma2(std::vector< std::tuple< ti_PDGID, td_Mass, TLorentzVector >> &v_prods)
Selection of the type of transformation matrix used in geometry.
TH2F * hist
Definition: plot.C:134
Int_t nbinsx
Definition: plot.C:23
std::string fPathToFile
CLHEP::HepRandomEngine & fEngine
#define MF_LOG_DEBUG(id)
std::vector< double > fT0
Beginning of time window to simulate in ns.
Representation of a node and its ancestry.
Definition: GeoNodePath.h:37
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.
Definition: GeometryCore.h:203
std::vector< double > fX0
Bottom corner x position (cm) in world coordinates.
TCEvent evt
Definition: DataStructs.cxx:8
Event generator information.
Definition: MCTruth.h:32
Float_t e
Definition: plot.C:35
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.
Definition: geo_vectors.h:229
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< double > fY0
Bottom corner y position (cm) in world coordinates.
void produce(art::Event &evt)