SingleGen_module.cc
Go to the documentation of this file.
1 ///
2 /// \brief Generator for specific kinematics - particularly single particles
3 /// \author brebel@fnal.gov
4 /// \date
5 ///
6 
7 // C++ includes.
8 #include <iostream>
9 #include <string>
10 #include <cmath>
11 #include <memory>
12 #include <vector>
13 
14 // ROOT includes
15 #include "TDatabasePDG.h"
16 #include "TFile.h"
17 #include "TH1F.h"
18 
19 // Framework includes
25 #include "fhiclcpp/ParameterSet.h"
30 
31 // NOvASoft includes
32 #include "Geometry/Geometry.h"
35 #include "nugen/EventGeneratorBase/evgenbase.h"
36 #include "SummaryData/RunData.h"
37 #include "SummaryData/SubRunData.h"
38 
39 #include "CLHEP/Random/RandFlat.h"
42 
43 
44 namespace evgen {
45  class SingleParticle;
46 
47  /// A module to check the results from the Monte Carlo generator
48  class SingleGen : public art::EDProducer {
49  public:
50  explicit SingleGen(fhicl::ParameterSet const& pset);
51  virtual ~SingleGen();
52 
53  void produce(art::Event& evt);
54  void reconfigure(const fhicl::ParameterSet& pset);
55  void beginRun(art::Run& run);
56  void endSubRun(art::SubRun& sr);
57 
58  private:
59 
60  void Sample(simb::MCTruth &truth);
61 
62  /// Depending on configuration of fPMeaning
63  /// calculate the momentum of the particle
64  double getMomentum(const double& pkf, const double& m) const;
65 
66  static const int kUNIF = 0;
67  static const int kGAUS = 1;
68  static const int kHIST = 2;
69 
70  std::vector<int> fPDG; ///< PDG code of particles to generate
71  int fPMeaning; ///< Meaning of P0, fSigmaP and fPDist.
72  ///< By default (fP0Meaning=0), P0 and sigmaP is momentum
73  ///< If fPMeaning=1, then P0 and sigmaP is kinetic energy in GeV
74  ///< If fPMeaning=2, then P0and sigmaP is total energy in GeV
75  int fCycle; ///< cycle number of MC production
76  std::vector<double> fP0; ///< Central momentum (GeV/c) to generate
77  std::vector<double> fSigmaP; ///< Variation in momenta (GeV/c)
78  std::vector<double> fSigmaP2; ///< Variation in momenta (GeV/c)
79  std::vector<int> fPDist; ///< How to distribute momenta (0=uniform, 1=gaussian)
80 
81  std::vector<double> fX0; ///< Central x position (cm)
82  std::vector<double> fY0; ///< Central y position (cm)
83  std::vector<double> fZ0; ///< Central z position (cm)
84  std::vector<double> fT0; ///< Central t position (ns)
85  std::vector<double> fSigmaX; ///< Variation in x position (cm)
86  std::vector<double> fSigmaY; ///< Variation in y position (cm)
87  std::vector<double> fSigmaZ; ///< Variation in z position (cm)
88  std::vector<double> fSigmaT; ///< Variation in t position (ns)
89  std::vector<double> fSigmaX2; ///< Variation in x position (cm)
90  std::vector<double> fSigmaY2; ///< Variation in y position (cm)
91  std::vector<double> fSigmaZ2; ///< Variation in z position (cm)
92  std::vector<double> fSigmaT2; ///< Variation in t position (ns)
93  std::vector<int> fPosDist; ///< How to distribute xyz (0=uniform, 1=gaussian)
94  std::vector<double> fCosZ0; ///< Cosine of central angle wrt z-axis
95  std::vector<double> fSigmaCosZ; ///< Size of variation of cosz
96  std::vector<double> fSigmaCosZ2; ///< Size of variation of cosz
97  std::vector<int> fCosZDist; ///< How to distribute cosz (0=uniform, 1=gaussian)
98  std::vector<double> fPhiXY0; ///< Central angle in the x-y plane (degrees)
99  std::vector<double> fSigmaPhiXY; ///< Size of variation in phixy (degrees)
100  std::vector<double> fSigmaPhiXY2; ///< Size of variation in phixy (degrees)
101  std::vector<int> fPhiDist; ///< How to distribute phixy (0=uniform, 1=gaussian)
102  std::vector<bool> fCoverage4pi; ///< Boolean to override other angle args and achieve 4pi coverage.
103 
105 
106  std::vector<float> sums;
107  float sum;
108  //int bins;
116 
117  TH1F *pmag_mu;
118  TH1F *theta_mu;
119  TH1F *phi_mu;
120  TH1F *pmag_pi;
121  TH1F *theta_pi;
122  TH1F *phi_pi;
123  TH1F *pmag_pr;
124  TH1F *theta_pr;
125  TH1F *phi_pr;
126  TH1F *pmag_ph;
127  TH1F *theta_ph;
128  TH1F *phi_ph;
129  TH1F *pmag_el;
130  TH1F *theta_el;
131  TH1F *phi_el;
132 
133  TH1F *posx_mu;
134  TH1F *posx_pi;
135  TH1F *posx_pr;
136  TH1F *posx_ph;
137  TH1F *posx_el;
138  TH1F *posy_mu;
139  TH1F *posy_pi;
140  TH1F *posy_pr;
141  TH1F *posy_ph;
142  TH1F *posy_el;
143  TH1F *posz_mu;
144  TH1F *posz_pi;
145  TH1F *posz_pr;
146  TH1F *posz_ph;
147  TH1F *posz_el;
148 
149  std::vector<float> sums_pmag_mu;
150  std::vector<float> sums_theta_mu;
151  std::vector<float> sums_phi_mu;
152  std::vector<float> sums_pmag_pi;
153  std::vector<float> sums_theta_pi;
154  std::vector<float> sums_phi_pi;
155  std::vector<float> sums_pmag_pr;
156  std::vector<float> sums_theta_pr;
157  std::vector<float> sums_phi_pr;
158  std::vector<float> sums_pmag_ph;
159  std::vector<float> sums_theta_ph;
160  std::vector<float> sums_phi_ph;
161  std::vector<float> sums_pmag_el;
162  std::vector<float> sums_theta_el;
163  std::vector<float> sums_phi_el;
164 
165  std::vector<float> sums_posx_mu;
166  std::vector<float> sums_posx_pi;
167  std::vector<float> sums_posx_pr;
168  std::vector<float> sums_posx_ph;
169  std::vector<float> sums_posx_el;
170  std::vector<float> sums_posy_mu;
171  std::vector<float> sums_posy_pi;
172  std::vector<float> sums_posy_pr;
173  std::vector<float> sums_posy_ph;
174  std::vector<float> sums_posy_el;
175  std::vector<float> sums_posz_mu;
176  std::vector<float> sums_posz_pi;
177  std::vector<float> sums_posz_pr;
178  std::vector<float> sums_posz_ph;
179  std::vector<float> sums_posz_el;
180 
196 
212 
213  // My TFS Hists
214  //TH1F *h_flatx, *h_flaty, *h_flatz;
215  //TH1F *h_pvecx, *h_pvecy, *h_pvecz;
216  };
217 }
218 
219 ////////////////////////////////////////////////////////////////////////
220 
221 namespace evgen{
222 
223  //____________________________________________________________________________
225  : fPMeaning(0)
226  , fCycle (pset.get<int>("Cycle", 0))
227  {
228  this->reconfigure(pset);
229 
230  // get the random number seed, use a random default if not specified
231  // in the configuration file.
232  int seed = pset.get< unsigned int >("Seed", evgb::GetRandomNumberSeed());
233 
234  createEngine( seed );
235 
236  produces< std::vector<simb::MCTruth> >();
237  produces< sumdata::SubRunData, art::InSubRun >();
238  produces< sumdata::RunData, art::InRun >();
239  }
240 
241  //____________________________________________________________________________
243 
244  //____________________________________________________________________________
246  {
247  fPDG = pset.get< std::vector<int> >("PDG");
248 
249  // Default to momentum
250  fPMeaning = pset.get< int >("PMeaning", 0);
251  if(fPMeaning == 1) {
252  std::cout<<"SingleGen: Using Kinetic energy for the meaning of P0, SigmaP and PDist\n";
253  }
254  else if(fPMeaning == 2) {
255  std::cout<<"SingleGen: Using Total energy for the meaning of P0, SigmaP and PDist\n";
256  }
257 
258  fP0 = pset.get< std::vector<double> >("P0");
259  fSigmaP = pset.get< std::vector<double> >("SigmaP");
260  fSigmaP2 = pset.get< std::vector<double> >("SigmaP2");
261  fPDist = pset.get< std::vector<int> >("PDist");
262  fX0 = pset.get< std::vector<double> >("X0");
263  fY0 = pset.get< std::vector<double> >("Y0");
264  fZ0 = pset.get< std::vector<double> >("Z0");
265  fT0 = pset.get< std::vector<double> >("T0");
266  fSigmaX = pset.get< std::vector<double> >("SigmaX");
267  fSigmaY = pset.get< std::vector<double> >("SigmaY");
268  fSigmaZ = pset.get< std::vector<double> >("SigmaZ");
269  fSigmaT = pset.get< std::vector<double> >("SigmaT");
270  fSigmaX2 = pset.get< std::vector<double> >("SigmaX2");
271  fSigmaY2 = pset.get< std::vector<double> >("SigmaY2");
272  fSigmaZ2 = pset.get< std::vector<double> >("SigmaZ2");
273  fSigmaT2 = pset.get< std::vector<double> >("SigmaT2");
274  fPosDist = pset.get< std::vector<int> >("PosDist");
275  fCosZ0 = pset.get< std::vector<double> >("CosZ0");
276  fSigmaCosZ = pset.get< std::vector<double> >("SigmaCosZ");
277  fSigmaCosZ2 = pset.get< std::vector<double> >("SigmaCosZ2");
278  fCosZDist = pset.get< std::vector<int> >("CosZDist");
279  fPhiXY0 = pset.get< std::vector<double> >("PhiXY0");
280  fSigmaPhiXY = pset.get< std::vector<double> >("SigmaPhiXY");
281  fSigmaPhiXY2= pset.get< std::vector<double> >("SigmaPhiXY2");
282  fPhiDist = pset.get< std::vector<int> >("PhiDist");
283  fCoverage4pi= pset.get< std::vector<bool> >("Coverage4pi");
284 
285  fROOTPath = pset.get< std::string >("ROOTPath");
286  }
287 
288  //____________________________________________________________________________
290  {
291  /*
292  art::ServiceHandle<art::TFileService> tfs;
293  h_flatx = tfs->make<TH1F>("UnitX", "Unit V Mom;X Number; Count" , 42, -1.1, 1.1);
294  h_flaty = tfs->make<TH1F>("UnitY", "Unit V Mom;Y Number; Count" , 42, -1.1, 1.1);
295  h_flatz = tfs->make<TH1F>("UnitZ", "Unit V Mom;Z Number; Count" , 42, -1.1, 1.1);
296 
297  h_pvecx = tfs->make<TH1F>("VecX", "Mom Vec;X Mom Vec; Count" , 100, -5, 5);
298  h_pvecy = tfs->make<TH1F>("VecY", "Mom Vec;Y Mom Vec; Count" , 100, -5, 5);
299  h_pvecz = tfs->make<TH1F>("VecZ", "Mom Vec;Z Mom Vec; Count" , 100, -5, 5);
300  */
301  // grab the geometry object to see what geometry we are using
303 
304  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetId(),
305  geo->FileBaseName(),
306  geo->ExtractGDML()));
307 
308  run.put(std::move(runcol));
309 
310  bool UseROOT = false;
311  for (unsigned int i=0; i<fPDG.size(); ++i)
312  if (fPDist[i] == kHIST || fPosDist[i] == kHIST || fCosZDist[i] == kHIST || fPhiDist[i] == kHIST)
313  UseROOT = true;
314 
315  if (UseROOT){
316  // Take histos from file
317  TFile f(fROOTPath.c_str());
318  assert(!f.IsZombie());
319  pmag_mu = (TH1F*)f.Get("pmag_mu");
320  theta_mu = (TH1F*)f.Get("theta_mu");
321  phi_mu = (TH1F*)f.Get("phi_mu");
322  pmag_pi = (TH1F*)f.Get("pmag_pi");
323  theta_pi = (TH1F*)f.Get("theta_pi");
324  phi_pi = (TH1F*)f.Get("phi_pi");
325  pmag_pr = (TH1F*)f.Get("pmag_pr");
326  theta_pr = (TH1F*)f.Get("theta_pr");
327  phi_pr = (TH1F*)f.Get("phi_pr");
328  pmag_ph = (TH1F*)f.Get("pmag_ph");
329  theta_ph = (TH1F*)f.Get("theta_ph");
330  phi_ph = (TH1F*)f.Get("phi_ph");
331  pmag_el = (TH1F*)f.Get("pmag_el");
332  theta_el = (TH1F*)f.Get("theta_el");
333  phi_el = (TH1F*)f.Get("phi_el");
334 
335  posx_mu = (TH1F*)f.Get("posx_mu");
336  posx_pi = (TH1F*)f.Get("posx_pi");
337  posx_pr = (TH1F*)f.Get("posx_pr");
338  posx_ph = (TH1F*)f.Get("posx_ph");
339  posx_el = (TH1F*)f.Get("posx_el");
340  posy_mu = (TH1F*)f.Get("posy_mu");
341  posy_pi = (TH1F*)f.Get("posy_pi");
342  posy_pr = (TH1F*)f.Get("posy_pr");
343  posy_ph = (TH1F*)f.Get("posy_ph");
344  posy_el = (TH1F*)f.Get("posy_el");
345  posz_mu = (TH1F*)f.Get("posz_mu");
346  posz_pi = (TH1F*)f.Get("posz_pi");
347  posz_pr = (TH1F*)f.Get("posz_pr");
348  posz_ph = (TH1F*)f.Get("posz_ph");
349  posz_el = (TH1F*)f.Get("posz_el");
350 
351  assert(pmag_mu);
352  assert(theta_mu);
353  assert(phi_mu);
354  assert(pmag_pi);
355  assert(theta_pi);
356  assert(phi_pi);
357  assert(pmag_pr);
358  assert(theta_pr);
359  assert(phi_pr);
360  assert(pmag_ph);
361  assert(theta_ph);
362  assert(phi_ph);
363  assert(pmag_el);
364  assert(theta_el);
365  assert(phi_el);
366 
367  assert(posx_mu);
368  assert(posx_pi);
369  assert(posx_pr);
370  assert(posx_ph);
371  assert(posx_el);
372  assert(posy_mu);
373  assert(posy_pi);
374  assert(posy_pr);
375  assert(posy_ph);
376  assert(posy_el);
377  assert(posz_mu);
378  assert(posz_pi);
379  assert(posz_pr);
380  assert(posz_ph);
381  assert(posz_el);
382 
383  bins_pmag_mu = pmag_mu->GetNbinsX();
384  bins_theta_mu = theta_mu->GetNbinsX();
385  bins_phi_mu = phi_mu->GetNbinsX();
386  bins_pmag_pi = pmag_pi->GetNbinsX();
387  bins_theta_pi = theta_pi->GetNbinsX();
388  bins_phi_pi = phi_pi->GetNbinsX();
389  bins_pmag_pr = pmag_pr->GetNbinsX();
390  bins_theta_pr = theta_pr->GetNbinsX();
391  bins_phi_pr = phi_pr->GetNbinsX();
392  bins_pmag_ph = pmag_ph->GetNbinsX();
393  bins_theta_ph = theta_ph->GetNbinsX();
394  bins_phi_ph = phi_ph->GetNbinsX();
395  bins_pmag_el = pmag_el->GetNbinsX();
396  bins_theta_el = theta_el->GetNbinsX();
397  bins_phi_el = phi_el->GetNbinsX();
398 
399  bins_posx_mu = posx_mu->GetNbinsX();
400  bins_posx_pi = posx_pi->GetNbinsX();
401  bins_posx_pr = posx_pr->GetNbinsX();
402  bins_posx_ph = posx_ph->GetNbinsX();
403  bins_posx_el = posx_el->GetNbinsX();
404  bins_posy_mu = posy_mu->GetNbinsX();
405  bins_posy_pi = posy_pi->GetNbinsX();
406  bins_posy_pr = posy_pr->GetNbinsX();
407  bins_posy_ph = posy_ph->GetNbinsX();
408  bins_posy_el = posy_el->GetNbinsX();
409  bins_posz_mu = posz_mu->GetNbinsX();
410  bins_posz_pi = posz_pi->GetNbinsX();
411  bins_posz_pr = posz_pr->GetNbinsX();
412  bins_posz_ph = posz_ph->GetNbinsX();
413  bins_posz_el = posz_el->GetNbinsX();
414 
415  // Take the CDFs
416  sum = 0;
417  for (Int_t i=0;i<bins_pmag_mu;i++){
418  sum = pmag_mu->GetBinContent(i) + sum;
419  sums_pmag_mu.push_back( sum * ( 1/pmag_mu->Integral() ) );
420  }
421  sum = 0;
422  for (Int_t i=0;i<bins_theta_mu;i++){
423  sum = theta_mu->GetBinContent(i) + sum;
424  sums_theta_mu.push_back( sum * ( 1/theta_mu->Integral() ) );
425  }
426  sum = 0;
427  for (Int_t i=0;i<bins_phi_mu;i++){
428  sum = phi_mu->GetBinContent(i) + sum;
429  sums_phi_mu.push_back( sum * ( 1/phi_mu->Integral() ) );
430  }
431 
432  sum = 0;
433  for (Int_t i=0;i<bins_pmag_pi;i++){
434  sum = pmag_pi->GetBinContent(i) + sum;
435  sums_pmag_pi.push_back( sum * ( 1/pmag_pi->Integral() ) );
436  }
437  sum = 0;
438  for (Int_t i=0;i<bins_theta_pi;i++){
439  sum = theta_pi->GetBinContent(i) + sum;
440  sums_theta_pi.push_back( sum * ( 1/theta_pi->Integral() ) );
441  }
442  sum = 0;
443  for (Int_t i=0;i<bins_phi_pi;i++){
444  sum = phi_pi->GetBinContent(i) + sum;
445  sums_phi_pi.push_back( sum * ( 1/phi_pi->Integral() ) );
446  }
447 
448  sum = 0;
449  for (Int_t i=0;i<bins_pmag_pr;i++){
450  sum = pmag_pr->GetBinContent(i) + sum;
451  sums_pmag_pr.push_back( sum * ( 1/pmag_pr->Integral() ) );
452  }
453  sum = 0;
454  for (Int_t i=0;i<bins_theta_pr;i++){
455  sum = theta_pr->GetBinContent(i) + sum;
456  sums_theta_pr.push_back( sum * ( 1/theta_pr->Integral() ) );
457  }
458  sum = 0;
459  for (Int_t i=0;i<bins_phi_pr;i++){
460  sum = phi_pr->GetBinContent(i) + sum;
461  sums_phi_pr.push_back( sum * ( 1/phi_pr->Integral() ) );
462  }
463 
464  sum = 0;
465  for (Int_t i=0;i<bins_pmag_ph;i++){
466  sum = pmag_ph->GetBinContent(i) + sum;
467  sums_pmag_ph.push_back( sum * ( 1/pmag_ph->Integral() ) );
468  }
469  sum = 0;
470  for (Int_t i=0;i<bins_theta_ph;i++){
471  sum = theta_ph->GetBinContent(i) + sum;
472  sums_theta_ph.push_back( sum * ( 1/theta_ph->Integral() ) );
473  }
474  sum = 0;
475  for (Int_t i=0;i<bins_phi_ph;i++){
476  sum = phi_ph->GetBinContent(i) + sum;
477  sums_phi_ph.push_back( sum * ( 1/phi_ph->Integral() ) );
478  }
479 
480  sum = 0;
481  for (Int_t i=0;i<bins_pmag_el;i++){
482  sum = pmag_el->GetBinContent(i) + sum;
483  sums_pmag_el.push_back( sum * ( 1/pmag_el->Integral() ) );
484  }
485  sum = 0;
486  for (Int_t i=0;i<bins_theta_el;i++){
487  sum = theta_el->GetBinContent(i) + sum;
488  sums_theta_el.push_back( sum * ( 1/theta_el->Integral() ) );
489  }
490  sum = 0;
491  for (Int_t i=0;i<bins_phi_el;i++){
492  sum = phi_el->GetBinContent(i) + sum;
493  sums_phi_el.push_back( sum * ( 1/phi_el->Integral() ) );
494  }
495 
496  sum = 0;
497  for (Int_t i=0;i<bins_posx_mu;i++){
498  sum = posx_mu->GetBinContent(i) + sum;
499  sums_posx_mu.push_back( sum * ( 1/posx_mu->Integral() ) );
500  }
501  sum = 0;
502  for (Int_t i=0;i<bins_posx_pi;i++){
503  sum = posx_pi->GetBinContent(i) + sum;
504  sums_posx_pi.push_back( sum * ( 1/posx_pi->Integral() ) );
505  }
506  sum = 0;
507  for (Int_t i=0;i<bins_posx_pr;i++){
508  sum = posx_pr->GetBinContent(i) + sum;
509  sums_posx_pr.push_back( sum * ( 1/posx_pr->Integral() ) );
510  }
511  sum = 0;
512  for (Int_t i=0;i<bins_posx_ph;i++){
513  sum = posx_ph->GetBinContent(i) + sum;
514  sums_posx_ph.push_back( sum * ( 1/posx_ph->Integral() ) );
515  }
516  sum = 0;
517  for (Int_t i=0;i<bins_posx_el;i++){
518  sum = posx_el->GetBinContent(i) + sum;
519  sums_posx_el.push_back( sum * ( 1/posx_el->Integral() ) );
520  }
521  sum = 0;
522  for (Int_t i=0;i<bins_posy_mu;i++){
523  sum = posy_mu->GetBinContent(i) + sum;
524  sums_posy_mu.push_back( sum * ( 1/posy_mu->Integral() ) );
525  }
526  sum = 0;
527  for (Int_t i=0;i<bins_posy_pi;i++){
528  sum = posy_pi->GetBinContent(i) + sum;
529  sums_posy_pi.push_back( sum * ( 1/posy_pi->Integral() ) );
530  }
531  sum = 0;
532  for (Int_t i=0;i<bins_posy_pr;i++){
533  sum = posy_pr->GetBinContent(i) + sum;
534  sums_posy_pr.push_back( sum * ( 1/posy_pr->Integral() ) );
535  }
536  sum = 0;
537  for (Int_t i=0;i<bins_posy_ph;i++){
538  sum = posy_ph->GetBinContent(i) + sum;
539  sums_posy_ph.push_back( sum * ( 1/posy_ph->Integral() ) );
540  }
541  sum = 0;
542  for (Int_t i=0;i<bins_posy_el;i++){
543  sum = posy_el->GetBinContent(i) + sum;
544  sums_posy_el.push_back( sum * ( 1/posy_el->Integral() ) );
545  }
546  sum = 0;
547  for (Int_t i=0;i<bins_posz_mu;i++){
548  sum = posz_mu->GetBinContent(i) + sum;
549  sums_posz_mu.push_back( sum * ( 1/posz_mu->Integral() ) );
550  }
551  sum = 0;
552  for (Int_t i=0;i<bins_posz_pi;i++){
553  sum = posz_pi->GetBinContent(i) + sum;
554  sums_posz_pi.push_back( sum * ( 1/posz_pi->Integral() ) );
555  }
556  sum = 0;
557  for (Int_t i=0;i<bins_posz_pr;i++){
558  sum = posz_pr->GetBinContent(i) + sum;
559  sums_posz_pr.push_back( sum * ( 1/posz_pr->Integral() ) );
560  }
561  sum = 0;
562  for (Int_t i=0;i<bins_posz_ph;i++){
563  sum = posz_ph->GetBinContent(i) + sum;
564  sums_posz_ph.push_back( sum * ( 1/posz_ph->Integral() ) );
565  }
566  sum = 0;
567  for (Int_t i=0;i<bins_posz_el;i++){
568  sum = posz_el->GetBinContent(i) + sum;
569  sums_posz_el.push_back( sum * ( 1/posz_el->Integral() ) );
570  }
571  }
572 
573  return;
574  }
575 
576  //____________________________________________________________________________
578  {
579  // store the cycle information
581  std::unique_ptr< sumdata::SubRunData > sd(new sumdata::SubRunData(fCycle));
582  sr.put(std::move(sd));
583 
584  }
585 
586  //____________________________________________________________________________
588  {
589 
590  ///unique_ptr allows ownership to be transferred to the art::Event after the ut statement
591  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
592 
593  simb::MCTruth truth;
595  Sample(truth);
596 
597  // std::cout << "put mctruth into the vector" << std::endl;
598  truthcol->push_back(truth);
599 
600  // std::cout << "add vector to the event " << truthcol->size() << std::endl;
601  evt.put(std::move(truthcol));
602 
603  return;
604  }
605 
606  //____________________________________________________________________________
608  {
609  // get the random number generator service and make some CLHEP generators
611  CLHEP::HepRandomEngine &engine = rng->getEngine();
612  CLHEP::RandFlat flat(engine);
613  CLHEP::RandGaussQ gauss(engine);
614 
615  /*
616  // Searchsorted for indices
617  //random_nums = gRandom->Uniform(0, 1);
618  random_nums = flat.fire(0, 1);
619  for (int j=1; j<bins_pmag_mu; j++){
620  if (random_nums<sums_pmag_mu.at(j) && random_nums>sums_pmag_mu.at(j-1))
621  random_index = j;
622  }
623  std::cout << "---------------------- HERE " << pmag_mu->GetBinCenter(random_index) << std::endl;
624  */
625 
626  //
627  // every event will have one of each particle species in the fPDG array
628  //
629  for (unsigned int i=0; i<fPDG.size(); ++i) {
630  // Momentum, kinetic energy or full energy
631  double pkf = 0.0;
632  if (fPDist[i] == kGAUS) {
633  pkf = abs( gauss.fire(fP0[i], fSigmaP[i]) );
634  }
635  else if (fPDist[i] == kUNIF) {
636  pkf = flat.fire(fSigmaP[i], fSigmaP2[i]);
637  }
638  else if (fPDist[i] == kHIST) {
639  random_nums_mag = flat.fire(0, 1);
640  switch (fPDG[i]) {
641  case 13: {
642  for (int j=1; j<bins_pmag_mu; j++)
644  random_index = j;
645  pkf = pmag_mu->GetBinCenter(random_index); } break;
646  case 211: {
647  for (int j=1; j<bins_pmag_pi; j++)
649  random_index = j;
650  pkf = pmag_pi->GetBinCenter(random_index); } break;
651  case 2212: {
652  for (int j=1; j<bins_pmag_pr; j++)
654  random_index = j;
655  pkf = pmag_pr->GetBinCenter(random_index); } break;
656  case 22: {
657  for (int j=1; j<bins_pmag_ph; j++)
659  random_index = j;
660  pkf = pmag_ph->GetBinCenter(random_index); } break;
661  case 11: {
662  for (int j=1; j<bins_pmag_el; j++)
664  random_index = j;
665  pkf = pmag_el->GetBinCenter(random_index); } break;
666  default : { std::cerr << "Not correct config" << std::endl; exit(-1); } break;
667  }
668  }
669  else {
670  std::cerr << "Not correct config: Mom" << std::endl;
671  exit(-1);
672  }
673 
674  // std::cout << "get the mass" << std::endl;
675  const TDatabasePDG* pdgt = TDatabasePDG::Instance();
676  const TParticlePDG* pdgp = pdgt->GetParticle(fPDG[i]);
677  // Mass in GeV
678  double m = 0.0;
679  if (pdgp) m = pdgp->Mass();
680 
681  // Momentum in GeV/c
682  const double p = getMomentum(pkf, m);
683 
684  // Choose position
685  double x[4];
686  if (fPosDist[i] == kGAUS) {
687  x[0] = gauss.fire(fX0[i], fSigmaX[i]);
688  x[1] = gauss.fire(fY0[i], fSigmaY[i]);
689  x[2] = gauss.fire(fZ0[i], fSigmaZ[i]);
690  x[3] = gauss.fire(fT0[i], fSigmaT[i]);
691  }
692  else if (fPosDist[i] == kUNIF) {
693  x[0] = flat.fire(fSigmaX[i], fSigmaX2[i]);
694  x[1] = flat.fire(fSigmaY[i], fSigmaY2[i]);
695  x[2] = flat.fire(fSigmaZ[i], fSigmaZ2[i]);
696  x[3] = flat.fire(fSigmaT[i], fSigmaT2[i]);
697  }
698  else if (fPosDist[i] == kHIST) {
699  x[0] = flat.fire(fSigmaX[i], fSigmaX2[i]);
700  x[1] = flat.fire(fSigmaY[i], fSigmaY2[i]);
701  x[2] = flat.fire(fSigmaZ[i], fSigmaZ2[i]);
702  x[3] = flat.fire(fSigmaT[i], fSigmaT2[i]);
703 
704  random_nums_posx = flat.fire(0, 1);
705  random_nums_posy = flat.fire(0, 1);
706  random_nums_posz = flat.fire(0, 1);
707  switch (fPDG[i]) {
708  case 13: {
709  for (int j=1; j<bins_posx_mu; j++)
711  random_index = j;
712  x[0] = posx_mu->GetBinCenter(random_index);
713  for (int j=1; j<bins_posy_mu; j++)
715  random_index = j;
716  x[1] = posy_mu->GetBinCenter(random_index);
717  for (int j=1; j<bins_posz_mu; j++)
719  random_index = j;
720  x[2] = posz_mu->GetBinCenter(random_index);
721  } break;
722  case 211: {
723  for (int j=1; j<bins_posx_pi; j++)
725  random_index = j;
726  x[0] = posx_pi->GetBinCenter(random_index);
727  for (int j=1; j<bins_posy_pi; j++)
729  random_index = j;
730  x[1] = posy_pi->GetBinCenter(random_index);
731  for (int j=1; j<bins_posz_pi; j++)
733  random_index = j;
734  x[2] = posz_pi->GetBinCenter(random_index);
735  } break;
736  case 2212: {
737  for (int j=1; j<bins_posx_pr; j++)
739  random_index = j;
740  x[0] = posx_pr->GetBinCenter(random_index);
741  for (int j=1; j<bins_posy_pr; j++)
743  random_index = j;
744  x[1] = posy_pr->GetBinCenter(random_index);
745  for (int j=1; j<bins_posz_pr; j++)
747  random_index = j;
748  x[2] = posz_pr->GetBinCenter(random_index);
749  } break;
750  case 22: {
751  for (int j=1; j<bins_posx_ph; j++)
753  random_index = j;
754  x[0] = posx_ph->GetBinCenter(random_index);
755  for (int j=1; j<bins_posy_ph; j++)
757  random_index = j;
758  x[1] = posy_ph->GetBinCenter(random_index);
759  for (int j=1; j<bins_posz_ph; j++)
761  random_index = j;
762  x[2] = posz_ph->GetBinCenter(random_index);
763  } break;
764  case 11: {
765  for (int j=1; j<bins_posx_el; j++)
767  random_index = j;
768  x[0] = posx_el->GetBinCenter(random_index);
769  for (int j=1; j<bins_posy_el; j++)
771  random_index = j;
772  x[1] = posy_el->GetBinCenter(random_index);
773  for (int j=1; j<bins_posz_el; j++)
775  random_index = j;
776  x[2] = posz_el->GetBinCenter(random_index);
777  } break;
778  default : { std::cerr << "Not correct config" << std::endl; exit(-1); } break;
779  }
780  }
781  else {
782  std::cerr << "Not correct config: Pos" << std::endl;
783  exit(-1);
784  }
785 
786  // std::cout << "set the position" << std::endl;
787  const TLorentzVector pos(x[0], x[1], x[2], x[3]);
788 
789  // Choose angles
790 
791  double cosz, phi;
792  unsigned int itry;
793  for (itry=0; itry<1000000; ++itry) {
794  if (fCosZDist[i] == kGAUS) {
795  cosz = gauss.fire(fCosZ0[i], fSigmaCosZ[i]);
796  }
797  else if (fCosZDist[i] == kUNIF) {
798  cosz = flat.fire(fSigmaCosZ[i], fSigmaCosZ2[i]);
799  }
800  else if (fCosZDist[i] == kHIST) {
801  //random_nums_theta = flat.fire(0, 1);
802  random_nums_theta = gRandom->Uniform(0, 1);
803  switch (fPDG[i]) {
804  case 13: {
805  for (int j=1; j<bins_theta_mu; j++)
807  random_index = j;
808  cosz = theta_mu->GetBinCenter(random_index);
809  } break;
810  case 211: {
811  for (int j=1; j<bins_theta_pi; j++)
813  random_index = j;
814  cosz = theta_pi->GetBinCenter(random_index);
815  } break;
816  case 2212: {
817  for (int j=1; j<bins_theta_pr; j++)
819  random_index = j;
820  cosz = theta_pr->GetBinCenter(random_index);
821  } break;
822  case 22: {
823  for (int j=1; j<bins_theta_ph; j++)
825  random_index = j;
826  cosz = theta_ph->GetBinCenter(random_index);
827  } break;
828  case 11: {
829  for (int j=1; j<bins_theta_el; j++)
831  random_index = j;
832  cosz = theta_el->GetBinCenter(random_index);
833  } break;
834  default : { std::cerr << "Not correct config" << std::endl; exit(-1); } break;
835  }
836  }
837  else {
838  std::cerr << "Not correct config: CosZ" << std::endl;
839  exit(-1);
840  }
841  if (cosz>=-1.0 && cosz<=1.0) break;
842  }
843 
844  if (cosz<-1.0 || cosz>1.0) {
845  mf::LogError("SingleGen") << __FILE__ << ":" << __LINE__
846  << " Failed to find an acceptable cos(theta_z)"
847  << " after many tries.\n"
848  << " Please adjust CosZ0 and SigmaCosZ0 in your"
849  << " SingleGen.fcl file and rerun";
850  abort();
851  }
852 
853  if (fPhiDist[i] == kGAUS) {
854  phi = gauss.fire(fPhiXY0[i], fSigmaPhiXY[i]);
855  }
856  else if(fPhiDist[i] == kUNIF) {
857  phi = flat.fire(fSigmaPhiXY[i], fSigmaPhiXY2[i]);
858  }
859  else if(fPhiDist[i] == kHIST) {
860  random_nums_phi = flat.fire(0, 1);
861  switch (fPDG[i]) {
862  case 13: {
863  for (int j=1; j<bins_phi_mu; j++)
865  random_index = j;
866  phi = phi_mu->GetBinCenter(random_index);
867  } break;
868  case 211: {
869  for (int j=1; j<bins_phi_pi; j++)
871  random_index = j;
872  phi = phi_pi->GetBinCenter(random_index);
873  } break;
874  case 2212: {
875  for (int j=1; j<bins_phi_pr; j++)
877  random_index = j;
878  phi = phi_pr->GetBinCenter(random_index);
879  } break;
880  case 22: {
881  for (int j=1; j<bins_phi_ph; j++)
883  random_index = j;
884  phi = phi_ph->GetBinCenter(random_index);
885  } break;
886  case 11: {
887  for (int j=1; j<bins_phi_el; j++)
889  random_index = j;
890  phi = phi_el->GetBinCenter(random_index);
891  } break;
892  default : { std::cerr << "Not correct config" << std::endl; exit(-1); } break;
893  }
894  }
895  else {
896  std::cerr << "Not correct config: Phi" << std::endl;
897  exit(-1);
898  }
899 
900 
901  const double sinz = sqrt(1.0-cosz*cosz);
902  const double sinphi = sin(phi*M_PI/180.0);
903  const double cosphi = cos(phi*M_PI/180.0);
904  //
905  // set track id to -i as these are all primary particles and have
906  // id <= 0
907  //
908  const int trackid = -1*(i+1);
909 
910  const std::string primary("primary");
911  simb::MCParticle part(trackid, fPDG[i], primary);
912 
913  //std::cout << "!!!!!!!!!!!!!! " << cosz << " " << p << std::endl;
914  TLorentzVector pvec;
915  if (fCoverage4pi[i]) {
916  // Make a unit vector, of gaussian sampled random numbers.
917  TVector3 flat_unit_p = TVector3(gauss.fire(), gauss.fire(), gauss.fire()).Unit();
918  pvec.SetPxPyPzE(p*flat_unit_p.Px(),
919  p*flat_unit_p.Py(),
920  p*flat_unit_p.Pz(),
921  sqrt(p*p+m*m));
922  /*
923  h_flatx -> Fill( flat_unit_p.Px() );
924  h_flaty -> Fill( flat_unit_p.Py() );
925  h_flatz -> Fill( flat_unit_p.Pz() );
926  h_pvecx -> Fill( pvec.Px() );
927  h_pvecy -> Fill( pvec.Py() );
928  h_pvecz -> Fill( pvec.Pz() );
929  */
930  } else {
931  pvec.SetPxPyPzE(cosphi*sinz*p,
932  sinphi*sinz*p,
933  cosz*p,
934  sqrt(p*p+m*m));
935  }
936 
937  part.AddTrajectoryPoint(pos, pvec);
938 
939  mct.Add(part);
940  } //end loop over particles
941  }
942 
943  //------------------------------------------------------------------------------
944  double SingleGen::getMomentum(const double& pkf, const double& m) const{
945 
946  if (fPMeaning == 0) return pkf;
947 
948  double total_energy = 0.0;
949  if (fPMeaning == 1){
950  total_energy = pkf + m;
951  }
952  else if(fPMeaning == 2){
953  total_energy = pkf;
954  }
955  else{
956  total_energy = sqrt(pkf*pkf + m*m);
957  }
958 
959  return sqrt(total_energy*total_energy - m*m);
960  }
961 
962  //------------------------------------------------------------------------------
964 
965 } //end namespace
966 
967 
std::vector< double > fSigmaT2
Variation in t position (ns)
std::vector< double > fSigmaP2
Variation in momenta (GeV/c)
std::vector< int > fCosZDist
How to distribute cosz (0=uniform, 1=gaussian)
std::vector< float > sums_phi_el
A module to check the results from the Monte Carlo generator.
std::vector< float > sums_posy_pi
std::vector< float > sums_theta_pi
std::vector< float > sums_posz_el
std::vector< float > sums_posz_pr
std::vector< double > fSigmaPhiXY2
Size of variation in phixy (degrees)
std::vector< float > sums
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:256
std::vector< float > sums_posx_mu
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
std::vector< float > sums_phi_pi
std::vector< double > fSigmaT
Variation in t position (ns)
std::vector< float > sums_posy_ph
double getMomentum(const double &pkf, const double &m) const
std::vector< double > fSigmaZ2
Variation in z position (cm)
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
std::vector< bool > fCoverage4pi
Boolean to override other angle args and achieve 4pi coverage.
base_engine_t & createEngine(seed_t seed)
std::vector< float > sums_posy_mu
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
std::vector< int > fPosDist
How to distribute xyz (0=uniform, 1=gaussian)
std::vector< double > fSigmaCosZ2
Size of variation of cosz.
std::vector< float > sums_theta_mu
std::vector< float > sums_posz_pi
std::vector< float > sums_phi_ph
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
static const int kUNIF
void reconfigure(const fhicl::ParameterSet &pset)
Particle class.
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
std::vector< float > sums_posx_el
std::vector< float > sums_phi_pr
#define M_PI
Definition: SbMath.h:34
std::vector< float > sums_theta_el
std::vector< float > sums_posz_ph
static const int kGAUS
std::vector< float > sums_posy_pr
Definition: Run.h:31
std::vector< float > sums_posx_pr
std::vector< double > fSigmaZ
Variation in z position (cm)
std::vector< float > sums_pmag_pi
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double sd(Eigen::VectorXd x)
std::vector< float > sums_pmag_ph
static const int kHIST
std::vector< double > fT0
Central t position (ns)
std::vector< double > fSigmaCosZ
Size of variation of cosz.
base_engine_t & getEngine() const
TString part[npart]
Definition: Style.C:32
std::vector< int > fPDist
How to distribute momenta (0=uniform, 1=gaussian)
std::vector< double > fSigmaPhiXY
Size of variation in phixy (degrees)
unsigned int seed
Definition: runWimpSim.h:102
std::vector< double > fP0
Central momentum (GeV/c) to generate.
std::vector< float > sums_theta_ph
single particles thrown at the detector
Definition: MCTruth.h:26
std::vector< float > sums_posx_ph
std::vector< float > sums_pmag_mu
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
std::vector< float > sums_pmag_el
std::vector< int > fPhiDist
How to distribute phixy (0=uniform, 1=gaussian)
static constexpr Double_t gauss
Definition: Munits.h:360
TVector3 Unit() const
std::vector< double > fCosZ0
Cosine of central angle wrt z-axis.
std::vector< float > sums_posx_pi
std::vector< double > fY0
Central y position (cm)
std::vector< int > fPDG
PDG code of particles to generate.
Definition: run.py:1
std::vector< double > fSigmaY2
Variation in y position (cm)
std::vector< float > sums_posz_mu
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::vector< double > fSigmaX2
Variation in x position (cm)
std::vector< double > fSigmaP
Variation in momenta (GeV/c)
void produce(art::Event &evt)
T sin(T number)
Definition: d0nt_math.hpp:132
std::vector< double > fZ0
Central z position (cm)
std::vector< double > fSigmaX
Variation in x position (cm)
int fCycle
cycle number of MC production
ProductID put(std::unique_ptr< PROD > &&)
void endSubRun(art::SubRun &sr)
void beginRun(art::Run &run)
exit(0)
std::vector< double > fX0
Central x position (cm)
std::vector< float > sums_pmag_pr
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
std::vector< double > fPhiXY0
Central angle in the x-y plane (degrees)
std::string ExtractGDML() const
Extract contents from fGDMLFile and return as a string.
void Sample(simb::MCTruth &truth)
Event generator information.
Definition: MCTruth.h:32
Helper for AttenCurve.
Definition: Path.h:10
Module to generate only pions from cosmic rays.
SingleGen(fhicl::ParameterSet const &pset)
Encapsulate the geometry of one entire detector (near, far, ndos)
std::vector< float > sums_phi_mu
std::string FileBaseName() const
std::vector< float > sums_theta_pr
std::vector< double > fSigmaY
Variation in y position (cm)
std::vector< float > sums_posy_el