MixSimEvents_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 // \file MixSimEvents_module.cc
3 // \author Jim Musser - jmusser@indiana.edu
4 // Gavin S. Davies - gsdavies@indiana.edu
5 // \brief Module for mixing MC files with MC files.
6 //
7 // Read zero or more events from an art event-data file and mix them
8 // into the current event. Part of the job is to update art::Ptr
9 // objects to point into the mixed collections. The simb::MCTruth
10 // simb::MCFlux objects are always read and mixed.
11 // NOvA-specific sim::Particle and sim::FLSHitList objects are also
12 // mixed.
13 //
14 // More detailed documentation on Product Mixing can be found here:
15 // https://cdcvs.fnal.gov/redmine/projects/art/wiki/Product_Mixing
16 //
17 ////////////////////////////////////////////////////////////////////////////
18 
19 // nusoft and NOVA includes
20 #include "nugen/EventGeneratorBase/evgenbase.h"
21 #include "Simulation/FLSHitList.h"
25 #include "Simulation/Particle.h"
26 #include "Simulation/TrueEnergy.h"
27 #include "SummaryData/SpillData.h"
28 #include "Utilities/AssociationUtil.h"
29 
30 #include "dk2nu/tree/dk2nu.h"
31 #include "dk2nu/tree/NuChoice.h"
32 
33 // Includes from art
42 
43 // ROOT includes
44 #include "TH1F.h"
45 
46 // Other third party includes
48 
49 // C++ includes
50 #include <memory>
51 #include <string>
52 #include <vector>
53 
54 
55 namespace mix
56 {
58  {
59  public:
60 
61  explicit MixSimEvents(const fhicl::ParameterSet& p,
62  art::MixHelper& helper);
63 
64 
65  // Set number of secondary events to read per primary event
66  size_t nSecondaries();
67 
68  // Mixing functions:
69  typedef std::vector<simb::MCTruth> MCTruthVector;
70  bool mixMCTruth(const std::vector<const MCTruthVector*>& in,
71  MCTruthVector& out,
72  const art::PtrRemapper& prm );
73 
74  typedef std::vector<simb::MCFlux> MCFluxVector;
75  bool mixMCFlux(const std::vector<const MCFluxVector*>& in,
76  MCFluxVector& out,
77  const art::PtrRemapper& prm);
78 
79  typedef std::vector<bsim::Dk2Nu> Dk2NuVector;
80  bool mixDk2Nu(const std::vector<const Dk2NuVector*>& in,
81  Dk2NuVector& out,
82  const art::PtrRemapper& prm);
83 
84  typedef std::vector<bsim::NuChoice> NuChoiceVector;
85  bool mixNuChoice(const std::vector<const NuChoiceVector*>& in,
86  NuChoiceVector& out,
87  const art::PtrRemapper& prm);
88 
89  typedef std::vector<simb::GTruth> GTruthVector;
90  bool mixGTruth(const std::vector<const GTruthVector*>& in,
91  GTruthVector& out,
92  const art::PtrRemapper& prm);
93 
94  typedef std::vector<sim::FLSHitList> FLSHitListVector;
95  bool mixFLSHitList(const std::vector<const FLSHitListVector*>& in,
96  FLSHitListVector& out,
97  const art::PtrRemapper& prm);
98 
99  typedef std::vector<sim::Particle> ParticleVector;
100  bool mixParticle(const std::vector<const ParticleVector*>& in,
101  ParticleVector& out,
102  const art::PtrRemapper& prm);
103 
104  typedef std::vector<sim::TrueEnergy> TrueEnergyVector;
105  bool mixTrueEnergy(const std::vector<const TrueEnergyVector*>& in,
106  TrueEnergyVector& out,
107  const art::PtrRemapper& prm);
108 
109  bool mixSpillData(const std::vector<const sumdata::SpillData*>& in,
111  const art::PtrRemapper& prm);
112 
113  private:
114 
115  // Run-time configurable members.
116  std::string fGeneratorLabel; /// Module label of the producer that made the generator-level information
117  std::string fGeant4Label; /// Module label of the producer that made the geantgen-level information
118 
119  bool fSecondCollectionIsCosmic; /// Is secondary simulated cosmic?
120  bool fSecondCollectionSpillData; /// Include secondary SpillData object?
121  bool fMixDk2Nu; /// Also mix in dk2nu/NuChoice objects? Default false
122 
123  std::string fMode; /// Mixing mode: useFixed, usePoisson, usePoissonMinus1, usePOT
124 
125  // useFixed
126  unsigned int fNFixed; /// Fixed number of events to mix in
127 
128  // usePoisson and usePoissonMinus1
129  double fMeanPoisson; /// The number of mix-in events to choose on each event, from poisson distribution mean
130 
131  // usePOT
132  double fPerPOT;
133  double fRefPOT;
134  double fPOT;
135 
136  // RNG
137  std::unique_ptr<CLHEP::RandPoissonQ> fRandPoisson;
138 
139  // Non-run-time configurable members.
140  TH1F* hNEvents; /// Histogram to record the distribution of generated events.
141  size_t fNMixActual; /// The actual number of mix-in events chosen on this event.
142 
143  // The offsets returned from flattening the GenParticle collections.
144  // These are needed in Ptr remapping operations.
145  std::vector<size_t> fGenOffsets; /// offset from flattening
146 
147  };
148 }
149 
150 // Some helper classes used only within this file.
151 namespace
152 {
153  void PoissonBins(double mean, double width,
154  int& nbins, double& xlow, double& xhigh)
155  {
156  if(mean < 0){
157  int n0 = static_cast<int>(std::floor(std::abs(mean)));
158  nbins = 3;
159  xlow = n0 - 1.;
160  xhigh = n0 + 1.;
161  }
162  else if(mean < 2){
163  nbins = 10;
164  xlow = 0.;
165  xhigh = 10.;
166  }
167  else if(mean < 25){
168  nbins = 50;
169  xlow = 0.;
170  xhigh = 50.;
171  }
172  else if(mean < 60){
173  nbins = 100;
174  xlow = 0.;
175  xhigh = 100.;
176  }
177  else{
178  int xl = static_cast<int>(std::floor(mean - width*std::sqrt(mean)));
179  int xh = static_cast<int>(std::ceil(mean + width*std::sqrt(mean)));
180  nbins = xh - xl;
181  xlow = xl;
182  xhigh = xh;
183  if(nbins > 100){
184  nbins = 100;
185  }
186  }
187  return;
188  }
189 } // end anonymous namespace
190 
191 
192 
193 namespace mix {
194 
196  art::MixHelper &helper):
197  // Run-time configurable parameters
198  fGeneratorLabel (pSet.get<std::string>("GeneratorLabel")),
199  fGeant4Label (pSet.get<std::string>("Geant4Label")),
200  fSecondCollectionIsCosmic (pSet.get<bool>("secondCollectionIsCosmic")),
201  fSecondCollectionSpillData(pSet.get<bool>("secondCollectionSpillData")),
202  fMixDk2Nu (pSet.get<bool>("MixDk2Nu")),
203  fMode (pSet.get<std::string>("mode")),
204  fNFixed (pSet.get<unsigned int>("fixed")),
205  fMeanPoisson (pSet.get<double>("meanPoisson")),
206  fPerPOT (pSet.get<double>("perPOT")),
207  fRefPOT (pSet.get<double>("refPOT")),
208  fPOT (pSet.get<double>("POT")),
209  fRandPoisson(nullptr),
210  // Non-run-time configurable
211  hNEvents(0),
212  fNMixActual(0)
213  {
215 
216  // Register MixOp operations; the callbacks are called in the order
217  // they were registered.
219  &MixSimEvents::mixMCTruth, *this );
220 
223  &MixSimEvents::mixMCFlux, *this );
224 
226  &MixSimEvents::mixGTruth, *this );
227  if(fMixDk2Nu){
229  &MixSimEvents::mixDk2Nu, *this);
231  &MixSimEvents::mixNuChoice, *this);
232  }
233 
234  }
236  &MixSimEvents::mixParticle, *this );
237 
239  &MixSimEvents::mixTrueEnergy, *this );
240 
242  &MixSimEvents::mixFLSHitList, *this );
243 
246  &MixSimEvents::mixSpillData, *this );
247  }
248 
250  int nbins;
251  double xlow, xhigh;
252 
253  if (fMode == "useFixed") {
254  nbins = 3;
255  xlow = fNFixed - 1;
256  xhigh = fNFixed + 1;
257  }
258  else {
259  // Determine the poisson mean based on mode
260  double thePoissonMean = 0;
261  if(fMode == "usePoisson"){
262  thePoissonMean = fMeanPoisson;
263  }
264  else if(fMode == "usePoissonMinus1"){
265  thePoissonMean = fMeanPoisson - 1;
266  }
267  else if (fMode == "usePOT") {
268  thePoissonMean = fPerPOT * fPOT / fRefPOT;
269  }
270  else {
271  std::cout << " Unknown mode " << fMode << " in overlay. Aborting" << std::endl;
272  std::abort();
273  }
274 
275  if (thePoissonMean < 0) {
276  std::cout << "Have a negative poisson mean of " << thePoissonMean << " when overlaying in mode " << fMode << std::endl;
277  std::abort();
278  }
279 
280  // Now actually throw the poisson random number
282  int dummy(0);
283  engine.setSeed(evgb::GetRandomNumberSeed(), dummy);
284  fRandPoisson = std::unique_ptr<CLHEP::RandPoissonQ>(new CLHEP::RandPoissonQ(engine, thePoissonMean));
285  PoissonBins(thePoissonMean, 4., nbins, xlow, xhigh);
286  }
287 
288  hNEvents = tfs->make<TH1F>("hNEvents", "Number of Mixed in Events",
289  nbins, xlow, xhigh);
290 
291  } // end mix::MixSimEvents::MixSimEvents
292 
293  //
294  //*********************MixSimEvents Methods****************
295  //
296 
297  //---------------------------------------------------------------------
299  {
300  fNMixActual = (fMode != "useFixed") ? fRandPoisson->fire() : fNFixed;
301  hNEvents->Fill(fNMixActual);
302  return fNMixActual;
303  }
304 
305  //---------------------------------------------------------------------
306  bool MixSimEvents::mixMCTruth(const std::vector<const MCTruthVector*>& in,
308  const art::PtrRemapper& prm)
309  {
310  // There are no Ptr's to update; just need to flatten.
312  return true;
313  }
314 
315  //---------------------------------------------------------------------
316  bool MixSimEvents::mixMCFlux(const std::vector<const MCFluxVector*>& in,
317  MCFluxVector& out,
318  const art::PtrRemapper& prm)
319  {
320  // There are no Ptr's to update; just need to flatten.
322  return true;
323  }
324 
325  //---------------------------------------------------------------------
326  bool MixSimEvents::mixDk2Nu(const std::vector<const Dk2NuVector*>& in,
327  Dk2NuVector& out,
328  const art::PtrRemapper& prm)
329  {
330  // There are no Ptr's to update; just need to flatten.
332  return true;
333  }
334 
335  //---------------------------------------------------------------------
336  bool MixSimEvents::mixNuChoice(const std::vector<const NuChoiceVector*>& in,
338  const art::PtrRemapper& prm)
339  {
340  // There are no Ptr's to update; just need to flatten.
342  return true;
343  }
344 
345  //---------------------------------------------------------------------
346  bool MixSimEvents::mixGTruth(const std::vector<const GTruthVector*>& in,
347  GTruthVector& out,
348  const art::PtrRemapper& prm)
349  {
350  // There are no Ptr's to update; just need to flatten.
352  return true;
353  }
354 
355  //---------------------------------------------------------------------
356  bool MixSimEvents::mixParticle(const std::vector<const ParticleVector*>& in,
358  const art::PtrRemapper& prm)
359  {
360  // There are no Ptr's to update; just need to flatten.
362  return true;
363  }
364 
365  //---------------------------------------------------------------------
366  bool MixSimEvents::mixTrueEnergy(const std::vector<const TrueEnergyVector*>& in,
368  const art::PtrRemapper& prm)
369  {
370  // There are no Ptr's to update; just need to flatten.
372  return true;
373  }
374 
375  //---------------------------------------------------------------------
376  bool MixSimEvents::mixFLSHitList(const std::vector<const FLSHitListVector*>& in,
378  const art::PtrRemapper& prm)
379  {
380  // There are no Ptr's to update; just need to flatten.
382  return true;
383  }
384 
385 
386  bool MixSimEvents::mixSpillData(const std::vector<const sumdata::SpillData*>& in,
388  const art::PtrRemapper& prm)
389  {
390  out = *in[0];
391  return true;
392  }
393 
395 
396 } // end namespace mix
397 ////////////////////////////////////////////////////////////////////////////
bool mixNuChoice(const std::vector< const NuChoiceVector * > &in, NuChoiceVector &out, const art::PtrRemapper &prm)
bool fMixDk2Nu
Include secondary SpillData object?
std::vector< bsim::NuChoice > NuChoiceVector
bool mixMCTruth(const std::vector< const MCTruthVector * > &in, MCTruthVector &out, const art::PtrRemapper &prm)
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
unsigned int GetRandomNumberSeed()
Definition: evgenbase.h:22
std::vector< simb::MCTruth > MCTruthVector
bool mixParticle(const std::vector< const ParticleVector * > &in, ParticleVector &out, const art::PtrRemapper &prm)
std::vector< bsim::Dk2Nu > Dk2NuVector
DEFINE_ART_MODULE(TestTMapFile)
bool mixDk2Nu(const std::vector< const Dk2NuVector * > &in, Dk2NuVector &out, const art::PtrRemapper &prm)
std::vector< sim::Particle > ParticleVector
float abs(float number)
Definition: d0nt_math.hpp:39
object containing MC flux information
std::vector< sim::TrueEnergy > TrueEnergyVector
bool mixGTruth(const std::vector< const GTruthVector * > &in, GTruthVector &out, const art::PtrRemapper &prm)
void flattenCollections(std::vector< COLLECTION const * > const &in, COLLECTION &out)
const int nbins
Definition: cellShifts.C:15
bool mixTrueEnergy(const std::vector< const TrueEnergyVector * > &in, TrueEnergyVector &out, const art::PtrRemapper &prm)
void declareMixOp(InputTag const &inputTag, MixFunc< PROD, OPROD > mixFunc, bool outputProduct=true)
Definition: MixHelper.h:428
bool mixMCFlux(const std::vector< const MCFluxVector * > &in, MCFluxVector &out, const art::PtrRemapper &prm)
std::vector< simb::MCFlux > MCFluxVector
size_t fNMixActual
Histogram to record the distribution of generated events.
std::vector< size_t > fGenOffsets
The actual number of mix-in events chosen on this event.
std::vector< simb::GTruth > GTruthVector
bool mixFLSHitList(const std::vector< const FLSHitListVector * > &in, FLSHitListVector &out, const art::PtrRemapper &prm)
OStream cout
Definition: OStream.cxx:6
unsigned int fNFixed
Mixing mode: useFixed, usePoisson, usePoissonMinus1, usePOT.
bool fSecondCollectionIsCosmic
Module label of the producer that made the geantgen-level information.
std::vector< sim::FLSHitList > FLSHitListVector
double fMeanPoisson
Fixed number of events to mix in.
T * make(ARGS...args) const
ifstream in
Definition: comparison.C:7
bool mixSpillData(const std::vector< const sumdata::SpillData * > &in, sumdata::SpillData &out, const art::PtrRemapper &prm)
std::unique_ptr< CLHEP::RandPoissonQ > fRandPoisson
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
bool fSecondCollectionSpillData
Is secondary simulated cosmic?
MixSimEvents(const fhicl::ParameterSet &p, art::MixHelper &helper)
std::string fMode
Also mix in dk2nu/NuChoice objects? Default false.
std::string fGeant4Label
Module label of the producer that made the generator-level information.
virtual void setSeed(long seed, int)=0
double fPerPOT
The number of mix-in events to choose on each event, from poisson distribution mean.
fvar< T > ceil(const fvar< T > &x)
Definition: ceil.hpp:11
enum BeamMode string