LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ActiveVolumeVertexSampler.cxx
Go to the documentation of this file.
1 
14 
16 
17 namespace {
18  constexpr int MAX_BOX_ITERATIONS = 10000;
19 }
20 
21 //------------------------------------------------------------------------------
24  rndm::NuRandomService& rand_service,
25  const geo::Geometry& geom,
26  const std::string& generator_name)
27  : fVertexType(vertex_type_t::kSampled), fGeneratorName(generator_name), fTPCDist(nullptr)
28 {
29  // Configure the algorithm using the FHiCL parameters
30  this->reconfigure(conf, geom);
31 
32  // Register the TPC sampling engine with the seed service. If you need the
33  // seed later, get it from the seed service using the value of the variable
34  // generator_name as the instance name.
35  rndm::NuRandomService::seed_t tpc_seed = rand_service.registerEngine(
36  [this](rndm::NuRandomService::EngineId const& /* unused */,
37  rndm::NuRandomService::seed_t lar_seed) -> void {
38  auto seed = static_cast<uint_fast64_t>(lar_seed);
39  // Use the obtained seed to prepare the random number engine. This is
40  // an attempt to do a decent job, but optimally accomplishing this can
41  // be tricky (see, for example,
42  // http://www.pcg-random.org/posts/cpp-seeding-surprises.html)
43  std::seed_seq seed_sequence{seed};
44  fTPCEngine.seed(seed_sequence);
45  },
47  conf.get_PSet(),
48  {"seed"});
49 
50  // TODO: resolve the other workaround mentioned in the MARLEYHelper
51  // class, then fix this as well
52  uint_fast64_t tpc_cast_seed = static_cast<uint_fast64_t>(tpc_seed);
53  std::seed_seq tpc_seed_sequence{tpc_cast_seed};
54  fTPCEngine.seed(tpc_seed_sequence);
55 }
56 
57 //------------------------------------------------------------------------------
59 {
60  // sample a new position if needed
62 
63  // Sample a TPC index using the active masses as weights
64  size_t tpc_index = (*fTPCDist)(fTPCEngine);
65 
66  // Get the dimensions of the chosen TPC's active volume
67  const auto& tpc = geom.TPC(geo::TPCID(0, tpc_index));
68  double minX = tpc.MinX();
69  double maxX = tpc.MaxX();
70  double minY = tpc.MinY();
71  double maxY = tpc.MaxY();
72  double minZ = tpc.MinZ();
73  double maxZ = tpc.MaxZ();
74  std::uniform_real_distribution<double>::param_type x_range(minX, maxX);
75  std::uniform_real_distribution<double>::param_type y_range(minY, maxY);
76  std::uniform_real_distribution<double>::param_type z_range(minZ, maxZ);
77 
78  // Sample a location uniformly over this volume
79  std::uniform_real_distribution<double> uniform_dist;
80  double x = uniform_dist(fTPCEngine, x_range);
81  double y = uniform_dist(fTPCEngine, y_range);
82  double z = uniform_dist(fTPCEngine, z_range);
83  MF_LOG_INFO("ActiveVolumeVertexSampler " + fGeneratorName)
84  << "Sampled primary vertex in TPC #" << tpc_index << ", x = " << x << ", y = " << y
85  << ", z = " << z;
86 
87  // Update the vertex position 4-vector
88  fVertexPosition.SetXYZT(x, y, z, 0.); // TODO: add time sampling
89  }
90  else if (fVertexType == vertex_type_t::kBox) {
91  bool ok = false;
92  int num_iterations = 0;
93 
94  // TODO: reduce code duplication here
95  // Get the range of coordinates covered by the box
96  std::uniform_real_distribution<double>::param_type x_range(fXmin, fXmax);
97  std::uniform_real_distribution<double>::param_type y_range(fYmin, fYmax);
98  std::uniform_real_distribution<double>::param_type z_range(fZmin, fZmax);
99 
100  // Sample a location uniformly over this volume
101  std::uniform_real_distribution<double> uniform_dist;
102 
103  double x = 0.;
104  double y = 0.;
105  double z = 0.;
106  while (!ok && num_iterations < MAX_BOX_ITERATIONS) {
107  x = uniform_dist(fTPCEngine, x_range);
108  y = uniform_dist(fTPCEngine, y_range);
109  z = uniform_dist(fTPCEngine, z_range);
110 
111  // If we've been asked to check whether the vertex is within the
112  // active volume of a TPC, verify that. Otherwise just accept the
113  // first vertex position that was sampled unconditionally.
114  if (fCheckActive) {
115  size_t num_tpcs = geom.NTPC();
116  for (auto const& tpc : geom.Iterate<geo::TPCGeo>(geo::CryostatID{0})) {
117  double minX = tpc.MinX();
118  double maxX = tpc.MaxX();
119  double minY = tpc.MinY();
120  double maxY = tpc.MaxY();
121  double minZ = tpc.MinZ();
122  double maxZ = tpc.MaxZ();
123  if (x >= minX && x <= maxX && y >= minY && y <= maxY && z >= minZ && z <= maxZ) {
124  ok = true;
125  break;
126  }
127  }
128  }
129  else
130  ok = true;
131  }
132 
133  if (!ok)
134  throw cet::exception("ActiveVolumeVertexSampler " + fGeneratorName)
135  << "Failed to sample a vertex within a TPC active volume after " << MAX_BOX_ITERATIONS
136  << " iterations";
137 
138  MF_LOG_INFO("ActiveVolumeVertexSampler " + fGeneratorName)
139  << "Sampled primary vertex at x = " << x << ", y = " << y << ", z = " << z;
140 
141  // Update the vertex position 4-vector
142  fVertexPosition.SetXYZT(x, y, z, 0.); // TODO: add time sampling
143  }
144 
145  // if we're using a fixed vertex position, we don't need to do any sampling
146  return fVertexPosition;
147 }
148 
149 //------------------------------------------------------------------------------
152  const geo::Geometry& geom)
153 {
154  auto type = conf().type_();
155  if (type == "sampled") {
156 
158 
159  // Get the active masses of each of the TPCs in the current geometry. Use
160  // them as weights for sampling a TPC to use for the primary vertex.
161  std::vector<double> tpc_masses;
162  size_t num_tpcs = geom.NTPC();
163  for (auto const& tpc : geom.Iterate<geo::TPCGeo>(geo::CryostatID{0})) {
164  // For each TPC, use its active mass (returned in kg) as its sampling weight
165  tpc_masses.push_back(tpc.ActiveMass());
166  }
167 
168  // Load the discrete distribution used to sample TPCs with the up-to-date
169  // set of weights
170  fTPCDist.reset(new std::discrete_distribution<size_t>(tpc_masses.begin(), tpc_masses.end()));
171  }
172  else if (type == "fixed") {
173 
175 
176  auto vertex_pos = conf().position_();
177  double Vx = vertex_pos.at(0);
178  double Vy = vertex_pos.at(1);
179  double Vz = vertex_pos.at(2);
180 
181  fVertexPosition.SetXYZT(Vx, Vy, Vz, 0.); // TODO: add time setting
182  }
183  else if (type == "box") {
184 
186  auto min_pos = conf().min_position_();
187  auto max_pos = conf().max_position_();
188 
189  fXmin = min_pos.at(0);
190  fYmin = min_pos.at(1);
191  fZmin = min_pos.at(2);
192 
193  fXmax = max_pos.at(0);
194  fYmax = max_pos.at(1);
195  fZmax = max_pos.at(2);
196 
197  // By default, don't check whether each point is in a TPC active volume
198  fCheckActive = false;
199  // If the user specified this optional parameter, use that instead
200  conf().check_active_(fCheckActive);
201  }
202 
203  else
204  throw cet::exception("ActiveVolumeVertexSampler " + fGeneratorName)
205  << "Invalid vertex type '" << type << "' requested. Allowed values are"
206  << " 'sampled' and 'fixed'";
207 }
Float_t x
Definition: compare.C:6
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
unsigned int NTPC(CryostatID const &cryoid=cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:686
double minY
Definition: plot_hist.C:9
Float_t y
Definition: compare.C:6
Double_t z
Definition: plot.C:276
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:90
Geometry information for a single TPC.
Definition: TPCGeo.h:36
ParameterSet const & get_PSet() const
Definition: Table.h:70
double maxY
Definition: plot_hist.C:10
TLorentzVector sample_vertex_pos(const geo::Geometry &geom)
seed_t registerEngine(SeedMaster_t::Seeder_t seeder, std::string const instance="", std::optional< seed_t > const seed=std::nullopt)
Registers an existing engine with rndm::NuRandomService.
ActiveVolumeVertexSampler(const fhicl::Table< Config > &conf, rndm::NuRandomService &rand_service, const geo::Geometry &geom, const std::string &generator_name)
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
Algorithm that samples vertex locations uniformly within the active volume of a detector. It is fully experiment-agnostic and multi-TPC aware.
long seed
Definition: chem4.cc:67
Identifier for a engine, made of module name and optional instance name.
Definition: EngineId.h:22
The geometry of one entire detector, as served by art.
Definition: Geometry.h:181
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
#define MF_LOG_INFO(category)
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
void reconfigure(const fhicl::Table< Config > &conf, const geo::Geometry &geom)
std::unique_ptr< std::discrete_distribution< size_t > > fTPCDist
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Encapsulate the construction of a single detector plane.
The data type to uniquely identify a cryostat.
Definition: geo_types.h:192
art::detail::EngineCreator::seed_t seed_t