SupernovaGen_module.cc
Go to the documentation of this file.
1 ///
2 /// \brief SNOwGLoBES-based Event Generator for Supernovae
3 /// \author Justin Vasel <jvasel@indiana.edu>
4 /// \author Matt Tamsett <m.tamsett@sussex.ac.uk>
5 ///
6 
7 // C++ includes
8 #include <cmath>
9 #include <iostream>
10 #include <memory>
11 #include <string>
12 #include <vector>
13 
14 // ROOT includes
15 #include "TDatabasePDG.h"
16 #include "TFile.h"
17 #include "TH2F.h"
18 #include "TRandom3.h"
19 
20 // CLHEP includes
21 #include "CLHEP/Random/RandFlat.h"
23 
24 // Framework includes
29 #include "fhiclcpp/ParameterSet.h"
34 
35 // NOvASoft includes
36 #include "nugen/EventGeneratorBase/evgenbase.h"
37 #include "Geometry/Geometry.h"
38 #include "NovaDAQConventions/DAQConventions.h"
41 #include "SummaryData/RunData.h"
42 #include "SummaryData/SubRunData.h"
43 
44 
45 namespace evgen {
46  class SupernovaGen : public art::EDProducer {
47  public:
48  explicit SupernovaGen(fhicl::ParameterSet const& p);
49  virtual ~SupernovaGen();
50 
51  void produce(art::Event& e);
52  void reconfigure(const fhicl::ParameterSet& p);
53  void beginRun(art::Run& r);
54  void endSubRun(art::SubRun& sr);
55 
56  private:
57 
58  void Generate(std::unique_ptr< std::vector<simb::MCTruth>> &truthcol);
59 
60  /// Depending on configuration of fPMeaning
61  /// calculate the momentum of the particle
62  double GetMomentum(const double& pkf, const double& m) const;
63 
64  static const int kUNIF = 0;
65  static const int kGAUS = 1;
66 
67  int fPDG; ///< PDG code of particles to generate
68  int fPMeaning; ///< Meaning of P0, fSigmaP and fPDist.
69  ///< By default (fP0Meaning=0), P0 and sigmaP is momentum
70  /// If fPMeaning=1, then P0 and sigmaP is kinetic energy in GeV
71  /// If fPMeaning=2, then P0and sigmaP is total energy in GeV
72  int fCycle; ///< MC production cycle
73  double fXlow; ///< Lower x position (cm)
74  double fXhigh; ///< Upper x position (cm)
75  double fYlow; ///< Lower y position (cm)
76  double fYhigh; ///< Upper y position (cm)
77  double fZlow; ///< Lower z position (cm)
78  double fZhigh; ///< Upper z position (cm)
79  double fTlow; ///< Lower t position (ns)
80  double fThigh; ///< Upper t position (ns)
81  double fCosZ0; ///< Cosine of central angle wrt z-axis
82  double fPhiXY0; ///< Central angle in the x-y plane (degrees)
83  double fSigmaCosZ; ///< Size of variation of cosz
84  double fSigmaPhiXY; ///< Size of variation in phixy (degrees)
85  int fAngleDist; ///< How to distribute angles (0=uniform, 1=gaussian)
86  std::string fPositronPDFFile; ///< ROOT file name where the positron PDFs are stored
87  bool fOnlyNoise; ///< Generate no signal
88  // Multiplicative scaling factor to change the distance of the SN
89  // the given rates are proportional to the distance^2 (10^2)
90  double fDistanceScaling; ///< Scaling factor to PDF values
91  art::EventNumber_t fEventPadding; ///< Number of empty events to place at beginning of file
92 
93  // Histograms
94  TH2F* fPositronPDF;
95 
97  double fMinPositronEnergy = 0.0007488; // GeV
98  double fMaxPositronEnergy = 0.09975; // GeV
99 
100  double fSigmaP;
101  int fTimeIndex = 0;
102  };
103 }
104 
105 ////////////////////////////////////////////////////////////////////////
106 
107 namespace evgen{
108 
109  //----------------------------------------------------------------------------
111  : fPMeaning(0)
112  , fCycle (p.get<int>("Cycle", 0))
113  {
114  this->reconfigure(p);
115 
116  // get the random number seed, use a random default if not specified
117  // in the configuration file.
118  int seed = p.get<unsigned int>("Seed", evgb::GetRandomNumberSeed());
119 
120  createEngine( seed );
121 
122  produces<std::vector<simb::MCTruth>>();
123  produces<sumdata::SubRunData, art::InSubRun>();
124  produces<sumdata::RunData, art::InRun >();
125  }
126 
127  //----------------------------------------------------------------------------
129 
130  //----------------------------------------------------------------------------
132  {
133  fPDG = p.get<int>("PDG");
134  try{
135  fPMeaning = p.get<int>("PMeaning");
136  if (fPMeaning == 1) {
137  std::cout << "SupernovaGen: Using Kinetic energy for the meaning of P0, SigmaP and PDist" << std::endl;
138  }
139  else if(fPMeaning == 2) {
140  std::cout << "SupernovaGen: Using Total energy for the meaning of P0, SigmaP and PDist" << std::endl;
141  }
142  }
143  catch(...){
144  fPMeaning = 0;
145  }
146 
147  fXlow = p.get<double> ("Xlow");
148  fYlow = p.get<double> ("Ylow");
149  fZlow = p.get<double> ("Zlow");
150  fXhigh = p.get<double> ("Xhigh");
151  fYhigh = p.get<double> ("Yhigh");
152  fZhigh = p.get<double> ("Zhigh");
153  fTlow = p.get<double> ("Tlow");
154  fThigh = p.get<double> ("Thigh");
155  fCosZ0 = p.get<double> ("CosZ0");
156  fPhiXY0 = p.get<double> ("PhiXY0");
157  fSigmaPhiXY = p.get<double> ("SigmaPhiXY");
158  fSigmaCosZ = p.get<double> ("SigmaCosZ");
159  fAngleDist = p.get<int> ("AngleDist");
160  fPositronPDFFile = p.get<std::string> ("PositronPDFFile");
161  fOnlyNoise = p.get<bool> ("OnlyNoise");
162  fDistanceScaling = p.get<double> ("DistanceScaling");
163 
164  fEventPadding = p.get<art::EventNumber_t>("EventPadding");
165 
166  TFile *root_file = new TFile(fPositronPDFFile.c_str());
167  root_file->GetObject("ibd", fPositronPDF);
168 
169  fNtimePoints = fPositronPDF->GetNbinsX();
170  fNenergyBins = fPositronPDF->GetNbinsY();
171 
173  }
174 
175  //----------------------------------------------------------------------------
177  {
178  // Get the detector information and push it into the event as RunData
180  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetId(),
181  geo->FileBaseName(),
182  geo->ExtractGDML()));
183  r.put(std::move(runcol));
184  return;
185  }
186 
187  //----------------------------------------------------------------------------
189  {
190  // store the cycle information
192  std::unique_ptr< sumdata::SubRunData > sd(new sumdata::SubRunData(fCycle));
193  sr.put(std::move(sd));
194  }
195 
196  //----------------------------------------------------------------------------
198  {
199  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
200 
201  // Only generate particles if we're done padding with empty events
202  if (e.event() > fEventPadding){
203  fTimeIndex++;
204 
205  // check that time index is within the available bounds
206  // TODO exit gracefully rather than crashing
207  if (fTimeIndex > fNtimePoints) return;
208 
209  Generate(truthcol);
210  }
211 
212  e.put(std::move(truthcol));
213  return;
214  }
215 
216  //----------------------------------------------------------------------------
217  void SupernovaGen::Generate(std::unique_ptr<std::vector<simb::MCTruth>> &truthcol)
218  {
219  // Set up the generators
221 
222  CLHEP::HepRandomEngine &engine = rng->getEngine();
223  CLHEP::RandFlat flat(engine);
224  CLHEP::RandGaussQ gauss(engine);
225 
226  //
227  // every event will have the positron distribution from the fPositronPDF array
228  //
229  int generatedPositrons = 0;
230 
231  for (int energyBin = 0; energyBin < fNenergyBins; ++energyBin) {
232  // the central momentum is given by the bin centre
233  double P0 = fMinPositronEnergy + (double(energyBin + 1) * 2. * fSigmaP);
234  double positronPDF = 0.;
235 
236  if (!fOnlyNoise)
237  positronPDF = fPositronPDF->GetBinContent(fTimeIndex, energyBin + 1);
238 
239  TRandom3 rnd(0); // zero seeds randomly
240 
241  int nPositrons = rnd.Poisson(positronPDF * fDistanceScaling);
242 
243  // Silence this output if no positrons produced
244  if (nPositrons > 0) {
245  std::cout << "\tEnergy bin[" << energyBin
246  << "] (" << P0
247  << " GeV), positron PDF: " << positronPDF
248  << std::endl;
249 
250  // Generate a Poisson-fluctated number of positrons based on this
251  std::cout << "\tGenerating " << nPositrons
252  << " positrons" << std::endl;
253  }
254 
255  for (int positronIdx = 0; positronIdx < nPositrons; ++positronIdx) {
256  generatedPositrons++;
257  simb::MCTruth mct;
259 
260  // Momentum, kinetic energy or full energy
261  double pkf = 0.0;
262  pkf = flat.fire(P0 - fSigmaP, P0 + fSigmaP);
263 
264  const TDatabasePDG* pdgt = TDatabasePDG::Instance();
265  const TParticlePDG* pdgp = pdgt->GetParticle(fPDG);
266 
267  // Mass in GeV
268  double m = 0.0;
269  if (pdgp) m = pdgp->Mass();
270 
271  // Momentum in GeV/c
272  const double p = GetMomentum(pkf, m);
273 
274  // Choose position
275  double x[4];
276  x[0] = flat.fire(fXlow, fXhigh);
277  x[1] = flat.fire(fYlow, fYhigh);
278  x[2] = flat.fire(fZlow, fZhigh);
279  x[3] = flat.fire(fTlow, fThigh);
280 
281  std::cout << "\t\tpositron[ " << generatedPositrons-1
282  << "]: p: " << p
283  << " GeV, x: " << x[0]
284  << " cm, y: " << x[1]
285  << " cm, z: " << x[2]
286  << " cm, t: " << x[3]
287  << " ns" << std::endl;
288 
289  const TLorentzVector pos(x[0], x[1], x[2], x[3]);
290 
291  // Choose angles
292  double cosz, phi;
293  unsigned int tryIdx;
294  for (tryIdx = 0; tryIdx < 1000000; ++tryIdx) {
295  if (fAngleDist == kGAUS) {
296  cosz = gauss.fire(fCosZ0, fSigmaCosZ);
297  phi = gauss.fire(fPhiXY0, fSigmaPhiXY);
298  } else {
299  cosz = flat.fire(fCosZ0 - fSigmaCosZ,
300  fCosZ0 + fSigmaCosZ);
301 
302  phi = flat.fire(fPhiXY0 - fSigmaPhiXY,
303  fPhiXY0 + fSigmaPhiXY);
304  }
305  if (cosz >= -1.0 && cosz <= 1.0) break;
306  }
307  if (cosz < -1.0 || cosz > 1.0) {
308  mf::LogError("SupernovaGen") << __FILE__ << ":" << __LINE__
309  << " Failed to find an acceptable cos(theta_z)"
310  << " after many tries.\n"
311  << " Please adjust CosZ0 and SigmaCosZ0 in your"
312  << " SupernovaGen.fcl file and rerun";
313  abort();
314  }
315 
316  const double sinz = sqrt(1.0 - cosz * cosz);
317  const double sinphi = sin(phi * M_PI / 180.0);
318  const double cosphi = cos(phi * M_PI / 180.0);
319 
320  // Set track id to -i as these are all primary particles and have id <= 0
321  const int trackid = -1 * (generatedPositrons + 1);
322 
323  const std::string primary("primary");
324  simb::MCParticle part(trackid, fPDG, primary);
325 
326  const TLorentzVector pvec(cosphi * sinz * p,
327  sinphi * sinz * p,
328  cosz * p,
329  sqrt(p * p + m * m));
330 
331  part.AddTrajectoryPoint(pos, pvec);
332 
333  mct.Add(part);
334  truthcol->push_back(mct);
335  } // end of loop over positrons
336  } // end loop over energy bins particles
337 
338  std::cout << "Generated a total of " << generatedPositrons
339  << " positrons for event " << fTimeIndex << std::endl;
340  }
341 
342  //------------------------------------------------------------------------------
343  double SupernovaGen::GetMomentum(const double& pkf, const double& m) const{
344 
345  if (fPMeaning == 0) return pkf;
346 
347  double total_energy = 0.0;
348  if (fPMeaning == 1){
349  total_energy = pkf + m;
350  }
351  else if(fPMeaning == 2){
352  total_energy = pkf;
353  }
354  else{
355  total_energy = sqrt(pkf * pkf + m * m);
356  }
357 
358  return sqrt(total_energy * total_energy - m * m);
359  }
360 
362 
363 } //end namespace
double fYlow
Lower y position (cm)
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:256
int fAngleDist
How to distribute angles (0=uniform, 1=gaussian)
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
static const int kUNIF
const char * p
Definition: xmltok.h:285
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
T sqrt(T number)
Definition: d0nt_math.hpp:156
unsigned int GetRandomNumberSeed()
Definition: evgenbase.h:22
base_engine_t & createEngine(seed_t seed)
double fYhigh
Upper y position (cm)
int fPDG
PDG code of particles to generate.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
#define M_PI
Definition: SbMath.h:34
Definition: Run.h:31
static const int kGAUS
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double fPhiXY0
Central angle in the x-y plane (degrees)
double sd(Eigen::VectorXd x)
void reconfigure(const fhicl::ParameterSet &p)
base_engine_t & getEngine() const
TString part[npart]
Definition: Style.C:32
double fSigmaCosZ
Size of variation of cosz.
unsigned int seed
Definition: runWimpSim.h:102
double fXhigh
Upper x position (cm)
single particles thrown at the detector
Definition: MCTruth.h:26
SupernovaGen(fhicl::ParameterSet const &p)
T get(std::string const &key) const
Definition: ParameterSet.h:231
void Generate(std::unique_ptr< std::vector< simb::MCTruth >> &truthcol)
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
double fZhigh
Upper z position (cm)
double fXlow
Lower x position (cm)
std::string fPositronPDFFile
ROOT file name where the positron PDFs are stored.
void endSubRun(art::SubRun &sr)
caf::StandardRecord * sr
EventNumber_t event() const
Definition: Event.h:67
static constexpr Double_t gauss
Definition: Munits.h:360
double GetMomentum(const double &pkf, const double &m) const
double fCosZ0
Cosine of central angle wrt z-axis.
OStream cout
Definition: OStream.cxx:6
double fDistanceScaling
Scaling factor to PDF values.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T sin(T number)
Definition: d0nt_math.hpp:132
bool fOnlyNoise
Generate no signal.
double fTlow
Lower t position (ns)
double fZlow
Lower z position (cm)
ProductID put(std::unique_ptr< PROD > &&)
double fSigmaPhiXY
Size of variation in phixy (degrees)
IDNumber_t< Level::Event > EventNumber_t
Definition: IDNumber.h:117
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
int fCycle
MC production cycle.
void produce(art::Event &e)
void beginRun(art::Run &r)
double fThigh
Upper t position (ns)
Event generator information.
Definition: MCTruth.h:32
Float_t e
Definition: plot.C:35
Helper for AttenCurve.
Definition: Path.h:10
Module to generate only pions from cosmic rays.
Encapsulate the geometry of one entire detector (near, far, ndos)
std::string FileBaseName() const
art::EventNumber_t fEventPadding
Number of empty events to place at beginning of file.