NueAr40CCGenerator.cxx
Go to the documentation of this file.
1 //=============================================================================
2 // NueAr40CCGenerator.cxx
3 //
4 // Gleb Sinev, Duke, 2015
5 //=============================================================================
6 
7 #include "NueAr40CCGenerator.h"
8 
9 // Framework includes
10 #include "cetlib/search_path.h"
11 #include "cetlib_except/exception.h"
12 #include "fhiclcpp/ParameterSet.h"
13 
14 // nusimdata includes
18 
19 // ROOT includes
20 #include "TMath.h"
21 #include "TFile.h"
22 #include "TGraph.h"
23 #include "TLorentzVector.h"
24 
25 #include "CLHEP/Random/RandFlat.h"
26 #include "CLHEP/Random/RandPoisson.h"
27 
28 // C++ includes
29 #include <cmath>
30 
31 namespace evgen {
32 
33  //----------------------------------------------------------------------------
34  // Constructor
36  const& parameterSet)
37  : fNumberOfLevels ( 73 )
38  , fNumberOfStartLevels ( 21 )
39  , fBranchingRatios (fNumberOfLevels )
40  , fDecayTo (fNumberOfLevels )
41  , fMonoenergeticNeutrinos
42  (parameterSet.get< bool >("MonoenergeticNeutrinos"))
43  , fNeutrinoEnergy
44  (parameterSet.get< double >("NeutrinoEnergy" ))
45  , fEnergySpectrumFileName
46  (parameterSet.get< std::string >("EnergySpectrumFileName"))
47  , fUsePoissonDistribution
48  (parameterSet.get< bool >("UsePoissonDistribution"))
49  , fAllowZeroNeutrinos
50  (parameterSet.get< bool >("AllowZeroNeutrinos" ))
51  , fNumberOfNeutrinos
52  (parameterSet.get< int >("NumberOfNeutrinos" ))
53  , fNeutrinoTimeBegin
54  (parameterSet.get< double >("NeutrinoTimeBegin" ))
55  , fNeutrinoTimeEnd
56  (parameterSet.get< double >("NeutrinoTimeEnd" ))
57  {
58 
59  fActiveVolume.push_back
60  (parameterSet.get< std::vector< double > >("ActiveVolume0"));
61  fActiveVolume.push_back
62  (parameterSet.get< std::vector< double > >("ActiveVolume1"));
63 
65 
67 
68  }
69 
70  //----------------------------------------------------------------------------
71  // Main routine
72  std::vector<simb::MCTruth> NueAr40CCGenerator::Generate(CLHEP::HepRandomEngine& engine)
73  {
74 
75  std::vector<simb::MCTruth> truths;
76 
77  int NumberOfNu = this->GetNumberOfNeutrinos(engine);
78  for(int i = 0; i < NumberOfNu; ++i) {
79 
80  simb::MCTruth truth;
82 
83  // Loop until at least one neutrino is simulated
84  while (!truth.NParticles()) {
85  CreateKinematicsVector(truth, engine);
86  }
87 
88  truths.push_back(truth);
89 
90  }
91 
92  return truths;
93 
94  }
95 
96  //----------------------------------------------------------------------------
97  // Return normalized direction cosines of isotropic vector
98  std::vector< double > NueAr40CCGenerator::GetIsotropicDirection
99  (CLHEP::HepRandomEngine& engine) const
100  {
101 
102  CLHEP::RandFlat randFlat(engine);
103 
104  std::vector< double > isotropicDirection;
105 
106  double phi = 2*TMath::Pi()*randFlat.fire();
107  double cosTheta = 2*randFlat.fire() - 1;
108  double theta = TMath::ACos(cosTheta);
109 
110  // x, y, z
111  isotropicDirection.push_back(cos(phi)*sin(theta));
112  isotropicDirection.push_back(sin(phi)*sin(theta));
113  isotropicDirection.push_back(cosTheta);
114 
115  return isotropicDirection;
116 
117  }
118 
119  //----------------------------------------------------------------------------
120  // Return a random vector with 3D coordinates inside the active LAr volume
121  std::vector< double > NueAr40CCGenerator::GetUniformPosition
122  (CLHEP::HepRandomEngine& engine) const
123  {
124 
125  CLHEP::RandFlat randFlat(engine);
126 
127  std::vector< double > position;
128 
129  position.push_back(randFlat.
130  fire(fActiveVolume.at(0).at(0), fActiveVolume.at(1).at(0)));
131  position.push_back(randFlat.
132  fire(fActiveVolume.at(0).at(1), fActiveVolume.at(1).at(1)));
133  position.push_back(randFlat.
134  fire(fActiveVolume.at(0).at(2), fActiveVolume.at(1).at(2)));
135 
136  return position;
137 
138  }
139 
140  //----------------------------------------------------------------------------
141  // Get number of neutrinos to generate
143  (CLHEP::HepRandomEngine& engine) const
144  {
145 
147  {
148  CLHEP::RandPoisson randPoisson(engine);
149  int N = randPoisson.fire(fNumberOfNeutrinos);
150  if(N == 0 && !fAllowZeroNeutrinos) N = 1;
151  return N;
152  }
153 
154  return fNumberOfNeutrinos;
155 
156  }
157 
158  //----------------------------------------------------------------------------
159  // Sample uniform distribution to get a neutrino interaction time
161  (CLHEP::HepRandomEngine& engine) const
162  {
163 
164  CLHEP::RandFlat randFlat(engine);
165 
166  return randFlat.fire(fNeutrinoTimeBegin, fNeutrinoTimeEnd);
167 
168  }
169 
170  //----------------------------------------------------------------------------
171  // Sample energy spectrum from fEnergyProbabilityMap
172  // or return a constant value
174  (CLHEP::HepRandomEngine& engine) const
175  {
176 
178 
179  CLHEP::RandFlat randFlat(engine);
180 
181  double neutrinoEnergy = 0.0;
182 
183  double randomNumber = randFlat.fire();
184 
185  // We need this to get a previous entry in the map
186  std::pair< double, double > previousPair;
187 
188  for (std::map< double, double >::const_iterator energyProbability =
189  fEnergyProbabilityMap.begin(); energyProbability !=
190  fEnergyProbabilityMap.end(); ++energyProbability)
191  {
192  if (randomNumber < energyProbability->second)
193  {
194  if (energyProbability != fEnergyProbabilityMap.begin())
195  {
196  neutrinoEnergy = energyProbability->first -
197  (energyProbability->second - randomNumber)*
198  (energyProbability->first - previousPair.first)/
199  (energyProbability->second - previousPair.second);
200  break;
201  }
202  else
203  {
204  neutrinoEnergy = energyProbability->first;
205  break;
206  }
207  }
208  previousPair = *energyProbability;
209  }
210 
211  return neutrinoEnergy;
212 
213  }
214 
215  //----------------------------------------------------------------------------
216  // Read a neutrino spectrum from a ROOT file
218  {
219 
220  cet::search_path searchPath("FW_SEARCH_PATH");
221  std::string directoryName = "SupernovaNeutrinos/" +
223 
224  std::string fullName;
225  searchPath.find_file(directoryName, fullName);
226 
227  if (fullName.empty())
228  throw cet::exception("NueAr40CCGenerator")
229  << "Cannot find file with neutrino energy spectrum "
230  << fullName << '\n';
231 
232  TFile energySpectrumFile(fullName.c_str(), "READ");
233 
234  std::string energySpectrumGraphName = "NueSpectrum";
235  TGraph *energySpectrumGraph =
236  dynamic_cast< TGraph* >(energySpectrumFile.Get
237  (energySpectrumGraphName.c_str()));
238 
239  double integral = 0.0;
240  int numberOfPoints = energySpectrumGraph->GetN();
241  double *energyValues = energySpectrumGraph->GetX();
242  double *fluxValues = energySpectrumGraph->GetY();
243  for (int point = 0; point < numberOfPoints; ++point)
244  integral += fluxValues[point];
245 
246  double probability = 0.0;
247  for (int point = 0; point < numberOfPoints; ++point)
248  {
249  probability += fluxValues[point]/integral;
250  fEnergyProbabilityMap.insert(std::make_pair(energyValues[point],
251  probability));
252  }
253 
254  }
255 
256  //----------------------------------------------------------------------------
257  // Fill vectors with quantities necessary to simulate gammas and electrons
259  {
260 
261  // The level data -- there's probably a better way to initialize
262  fBranchingRatios.at(0).push_back(1.00);
263  fDecayTo.at(0).push_back(1);
264 
265  fBranchingRatios.at(1).push_back(1.00);
266  fDecayTo.at(1).push_back(4);
267 
268  fBranchingRatios.at(2).push_back(0.133891);
269  fBranchingRatios.at(2).push_back(0.41841);
270  fBranchingRatios.at(2).push_back(0.351464);
271  fBranchingRatios.at(2).push_back(0.0962343);
272  fDecayTo.at(2).push_back(38);
273  fDecayTo.at(2).push_back(43);
274  fDecayTo.at(2).push_back(59);
275  fDecayTo.at(2).push_back(72);
276 
277  fBranchingRatios.at(3).push_back(0.391705);
278  fBranchingRatios.at(3).push_back(0.460829);
279  fBranchingRatios.at(3).push_back(0.147465);
280  fDecayTo.at(3).push_back(61);
281  fDecayTo.at(3).push_back(65);
282  fDecayTo.at(3).push_back(71);
283 
284  fBranchingRatios.at(4).push_back(0.358974);
285  fBranchingRatios.at(4).push_back(0.641026);
286  fDecayTo.at(4).push_back(13);
287  fDecayTo.at(4).push_back(58);
288 
289  fBranchingRatios.at(5).push_back(0.247808);
290  fBranchingRatios.at(5).push_back(0.0468929);
291  fBranchingRatios.at(5).push_back(0.213496);
292  fBranchingRatios.at(5).push_back(0.11056);
293  fBranchingRatios.at(5).push_back(0.381243);
294  fDecayTo.at(5).push_back(40);
295  fDecayTo.at(5).push_back(52);
296  fDecayTo.at(5).push_back(67);
297  fDecayTo.at(5).push_back(71);
298  fDecayTo.at(5).push_back(72);
299 
300  fBranchingRatios.at(6).push_back(0.361025);
301  fBranchingRatios.at(6).push_back(0.056677);
302  fBranchingRatios.at(6).push_back(0.194099);
303  fBranchingRatios.at(6).push_back(0.388199);
304  fDecayTo.at(6).push_back(52);
305  fDecayTo.at(6).push_back(55);
306  fDecayTo.at(6).push_back(63);
307  fDecayTo.at(6).push_back(68);
308 
309  fBranchingRatios.at(7).push_back(0.0300963);
310  fBranchingRatios.at(7).push_back(0.0613965);
311  fBranchingRatios.at(7).push_back(0.0152488);
312  fBranchingRatios.at(7).push_back(0.0321027);
313  fBranchingRatios.at(7).push_back(0.0513644);
314  fBranchingRatios.at(7).push_back(0.0682183);
315  fBranchingRatios.at(7).push_back(0.073435);
316  fBranchingRatios.at(7).push_back(0.0100321);
317  fBranchingRatios.at(7).push_back(0.0120385);
318  fBranchingRatios.at(7).push_back(0.0922953);
319  fBranchingRatios.at(7).push_back(0.152488);
320  fBranchingRatios.at(7).push_back(0.401284);
321  fDecayTo.at(7).push_back(17);
322  fDecayTo.at(7).push_back(25);
323  fDecayTo.at(7).push_back(27);
324  fDecayTo.at(7).push_back(32);
325  fDecayTo.at(7).push_back(49);
326  fDecayTo.at(7).push_back(54);
327  fDecayTo.at(7).push_back(56);
328  fDecayTo.at(7).push_back(62);
329  fDecayTo.at(7).push_back(63);
330  fDecayTo.at(7).push_back(67);
331  fDecayTo.at(7).push_back(68);
332  fDecayTo.at(7).push_back(70);
333 
334  fBranchingRatios.at(8).push_back(0.0089877);
335  fBranchingRatios.at(8).push_back(0.06386);
336  fBranchingRatios.at(8).push_back(0.13245);
337  fBranchingRatios.at(8).push_back(0.473037);
338  fBranchingRatios.at(8).push_back(0.321665);
339  fDecayTo.at(8).push_back(43);
340  fDecayTo.at(8).push_back(56);
341  fDecayTo.at(8).push_back(67);
342  fDecayTo.at(8).push_back(70);
343  fDecayTo.at(8).push_back(71);
344 
345  fBranchingRatios.at(9).push_back(0.16129);
346  fBranchingRatios.at(9).push_back(0.193548);
347  fBranchingRatios.at(9).push_back(0.645161);
348  fDecayTo.at(9).push_back(37);
349  fDecayTo.at(9).push_back(65);
350  fDecayTo.at(9).push_back(72);
351 
352  fBranchingRatios.at(10).push_back(0.245826);
353  fBranchingRatios.at(10).push_back(0.0816327);
354  fBranchingRatios.at(10).push_back(0.20872);
355  fBranchingRatios.at(10).push_back(0.463822);
356  fDecayTo.at(10).push_back(40);
357  fDecayTo.at(10).push_back(56);
358  fDecayTo.at(10).push_back(60);
359  fDecayTo.at(10).push_back(71);
360 
361  fBranchingRatios.at(11).push_back(0.170984);
362  fBranchingRatios.at(11).push_back(0.310881);
363  fBranchingRatios.at(11).push_back(0.518135);
364  fDecayTo.at(11).push_back(29);
365  fDecayTo.at(11).push_back(54);
366  fDecayTo.at(11).push_back(66);
367 
368  fBranchingRatios.at(12).push_back(0.242424);
369  fBranchingRatios.at(12).push_back(0.757576);
370  fDecayTo.at(12).push_back(54);
371  fDecayTo.at(12).push_back(62);
372 
373  fBranchingRatios.at(13).push_back(0.159664);
374  fBranchingRatios.at(13).push_back(0.840336);
375  fDecayTo.at(13).push_back(48);
376  fDecayTo.at(13).push_back(58);
377 
378  fBranchingRatios.at(14).push_back(0.288136);
379  fBranchingRatios.at(14).push_back(0.145763);
380  fBranchingRatios.at(14).push_back(0.118644);
381  fBranchingRatios.at(14).push_back(0.108475);
382  fBranchingRatios.at(14).push_back(0.338983);
383  fDecayTo.at(14).push_back(56);
384  fDecayTo.at(14).push_back(66);
385  fDecayTo.at(14).push_back(70);
386  fDecayTo.at(14).push_back(71);
387  fDecayTo.at(14).push_back(72);
388 
389  fBranchingRatios.at(15).push_back(0.00869263);
390  fBranchingRatios.at(15).push_back(0.087274);
391  fBranchingRatios.at(15).push_back(0.0973574);
392  fBranchingRatios.at(15).push_back(0.288595);
393  fBranchingRatios.at(15).push_back(0.347705);
394  fBranchingRatios.at(15).push_back(0.170376);
395  fDecayTo.at(15).push_back(40);
396  fDecayTo.at(15).push_back(64);
397  fDecayTo.at(15).push_back(65);
398  fDecayTo.at(15).push_back(68);
399  fDecayTo.at(15).push_back(70);
400  fDecayTo.at(15).push_back(71);
401 
402  fBranchingRatios.at(16).push_back(1);
403  fDecayTo.at(16).push_back(65);
404 
405  fBranchingRatios.at(17).push_back(0.341763);
406  fBranchingRatios.at(17).push_back(0.0567327);
407  fBranchingRatios.at(17).push_back(0.164046);
408  fBranchingRatios.at(17).push_back(0.239234);
409  fBranchingRatios.at(17).push_back(0.198223);
410  fDecayTo.at(17).push_back(35);
411  fDecayTo.at(17).push_back(39);
412  fDecayTo.at(17).push_back(51);
413  fDecayTo.at(17).push_back(67);
414  fDecayTo.at(17).push_back(69);
415 
416  fBranchingRatios.at(18).push_back(0.106406);
417  fBranchingRatios.at(18).push_back(0.254103);
418  fBranchingRatios.at(18).push_back(0.0465855);
419  fBranchingRatios.at(18).push_back(0.529381);
420  fBranchingRatios.at(18).push_back(0.0635257);
421  fDecayTo.at(18).push_back(60);
422  fDecayTo.at(18).push_back(61);
423  fDecayTo.at(18).push_back(63);
424  fDecayTo.at(18).push_back(70);
425  fDecayTo.at(18).push_back(72);
426 
427  fBranchingRatios.at(19).push_back(0.0905797);
428  fBranchingRatios.at(19).push_back(0.228261);
429  fBranchingRatios.at(19).push_back(0.181159);
430  fBranchingRatios.at(19).push_back(0.137681);
431  fBranchingRatios.at(19).push_back(0.362319);
432  fDecayTo.at(19).push_back(43);
433  fDecayTo.at(19).push_back(45);
434  fDecayTo.at(19).push_back(52);
435  fDecayTo.at(19).push_back(70);
436  fDecayTo.at(19).push_back(71);
437 
438  fBranchingRatios.at(20).push_back(0.0316414);
439  fBranchingRatios.at(20).push_back(0.0415293);
440  fBranchingRatios.at(20).push_back(0.0362558);
441  fBranchingRatios.at(20).push_back(0.0481213);
442  fBranchingRatios.at(20).push_back(0.0903098);
443  fBranchingRatios.at(20).push_back(0.0929466);
444  fBranchingRatios.at(20).push_back(0.659196);
445  fDecayTo.at(20).push_back(30);
446  fDecayTo.at(20).push_back(32);
447  fDecayTo.at(20).push_back(46);
448  fDecayTo.at(20).push_back(61);
449  fDecayTo.at(20).push_back(64);
450  fDecayTo.at(20).push_back(66);
451  fDecayTo.at(20).push_back(70);
452 
453  fBranchingRatios.at(21).push_back(0.310757);
454  fBranchingRatios.at(21).push_back(0.398406);
455  fBranchingRatios.at(21).push_back(0.290837);
456  fDecayTo.at(21).push_back(64);
457  fDecayTo.at(21).push_back(66);
458  fDecayTo.at(21).push_back(70);
459 
460  fBranchingRatios.at(22).push_back(0.152542);
461  fBranchingRatios.at(22).push_back(0.847458);
462  fDecayTo.at(22).push_back(67);
463  fDecayTo.at(22).push_back(71);
464 
465  fBranchingRatios.at(23).push_back(0.168367);
466  fBranchingRatios.at(23).push_back(0.321429);
467  fBranchingRatios.at(23).push_back(0.510204);
468  fDecayTo.at(23).push_back(52);
469  fDecayTo.at(23).push_back(70);
470  fDecayTo.at(23).push_back(71);
471 
472  fBranchingRatios.at(24).push_back(0.0276437);
473  fBranchingRatios.at(24).push_back(0.078982);
474  fBranchingRatios.at(24).push_back(0.0245722);
475  fBranchingRatios.at(24).push_back(0.162352);
476  fBranchingRatios.at(24).push_back(0.184291);
477  fBranchingRatios.at(24).push_back(0.438789);
478  fBranchingRatios.at(24).push_back(0.0833699);
479  fDecayTo.at(24).push_back(36);
480  fDecayTo.at(24).push_back(53);
481  fDecayTo.at(24).push_back(62);
482  fDecayTo.at(24).push_back(64);
483  fDecayTo.at(24).push_back(70);
484  fDecayTo.at(24).push_back(71);
485  fDecayTo.at(24).push_back(72);
486 
487  fBranchingRatios.at(25).push_back(0.0163934);
488  fBranchingRatios.at(25).push_back(0.336276);
489  fBranchingRatios.at(25).push_back(0.226986);
490  fBranchingRatios.at(25).push_back(0.420345);
491  fDecayTo.at(25).push_back(43);
492  fDecayTo.at(25).push_back(67);
493  fDecayTo.at(25).push_back(68);
494  fDecayTo.at(25).push_back(70);
495 
496  fBranchingRatios.at(26).push_back(0.184332);
497  fBranchingRatios.at(26).push_back(0.0460829);
498  fBranchingRatios.at(26).push_back(0.460829);
499  fBranchingRatios.at(26).push_back(0.078341);
500  fBranchingRatios.at(26).push_back(0.230415);
501  fDecayTo.at(26).push_back(53);
502  fDecayTo.at(26).push_back(54);
503  fDecayTo.at(26).push_back(60);
504  fDecayTo.at(26).push_back(61);
505  fDecayTo.at(26).push_back(71);
506 
507  fBranchingRatios.at(27).push_back(0.0147145);
508  fBranchingRatios.at(27).push_back(0.0170689);
509  fBranchingRatios.at(27).push_back(0.0500294);
510  fBranchingRatios.at(27).push_back(0.329606);
511  fBranchingRatios.at(27).push_back(0.588582);
512  fDecayTo.at(27).push_back(36);
513  fDecayTo.at(27).push_back(46);
514  fDecayTo.at(27).push_back(56);
515  fDecayTo.at(27).push_back(67);
516  fDecayTo.at(27).push_back(68);
517 
518  fBranchingRatios.at(28).push_back(1);
519  fDecayTo.at(28).push_back(70);
520 
521  fBranchingRatios.at(29).push_back(0.280702);
522  fBranchingRatios.at(29).push_back(0.0935673);
523  fBranchingRatios.at(29).push_back(0.0409357);
524  fBranchingRatios.at(29).push_back(0.584795);
525  fDecayTo.at(29).push_back(63);
526  fDecayTo.at(29).push_back(66);
527  fDecayTo.at(29).push_back(68);
528  fDecayTo.at(29).push_back(70);
529 
530  fBranchingRatios.at(30).push_back(0.393701);
531  fBranchingRatios.at(30).push_back(0.283465);
532  fBranchingRatios.at(30).push_back(0.188976);
533  fBranchingRatios.at(30).push_back(0.133858);
534  fDecayTo.at(30).push_back(61);
535  fDecayTo.at(30).push_back(67);
536  fDecayTo.at(30).push_back(71);
537  fDecayTo.at(30).push_back(72);
538 
539  fBranchingRatios.at(31).push_back(0.0477737);
540  fBranchingRatios.at(31).push_back(0.194805);
541  fBranchingRatios.at(31).push_back(0.0245826);
542  fBranchingRatios.at(31).push_back(0.269017);
543  fBranchingRatios.at(31).push_back(0.463822);
544  fDecayTo.at(31).push_back(45);
545  fDecayTo.at(31).push_back(60);
546  fDecayTo.at(31).push_back(65);
547  fDecayTo.at(31).push_back(71);
548  fDecayTo.at(31).push_back(72);
549 
550  fBranchingRatios.at(32).push_back(0.105769);
551  fBranchingRatios.at(32).push_back(0.129808);
552  fBranchingRatios.at(32).push_back(0.0528846);
553  fBranchingRatios.at(32).push_back(0.480769);
554  fBranchingRatios.at(32).push_back(0.230769);
555  fDecayTo.at(32).push_back(46);
556  fDecayTo.at(32).push_back(56);
557  fDecayTo.at(32).push_back(66);
558  fDecayTo.at(32).push_back(70);
559  fDecayTo.at(32).push_back(71);
560 
561  fBranchingRatios.at(33).push_back(0.21875);
562  fBranchingRatios.at(33).push_back(0.78125);
563  fDecayTo.at(33).push_back(67);
564  fDecayTo.at(33).push_back(71);
565 
566  fBranchingRatios.at(34).push_back(0.102011);
567  fBranchingRatios.at(34).push_back(0.179598);
568  fBranchingRatios.at(34).push_back(0.718391);
569  fDecayTo.at(34).push_back(43);
570  fDecayTo.at(34).push_back(61);
571  fDecayTo.at(34).push_back(66);
572 
573  fBranchingRatios.at(35).push_back(0.00826446);
574  fBranchingRatios.at(35).push_back(0.393546);
575  fBranchingRatios.at(35).push_back(0.334514);
576  fBranchingRatios.at(35).push_back(0.263676);
577  fDecayTo.at(35).push_back(64);
578  fDecayTo.at(35).push_back(67);
579  fDecayTo.at(35).push_back(68);
580  fDecayTo.at(35).push_back(70);
581 
582  fBranchingRatios.at(36).push_back(0.056338);
583  fBranchingRatios.at(36).push_back(0.704225);
584  fBranchingRatios.at(36).push_back(0.239437);
585  fDecayTo.at(36).push_back(51);
586  fDecayTo.at(36).push_back(70);
587  fDecayTo.at(36).push_back(71);
588 
589  fBranchingRatios.at(37).push_back(0.21875);
590  fBranchingRatios.at(37).push_back(0.78125);
591  fDecayTo.at(37).push_back(67);
592  fDecayTo.at(37).push_back(70);
593 
594  fBranchingRatios.at(38).push_back(0.181818);
595  fBranchingRatios.at(38).push_back(0.757576);
596  fBranchingRatios.at(38).push_back(0.0606061);
597  fDecayTo.at(38).push_back(66);
598  fDecayTo.at(38).push_back(71);
599  fDecayTo.at(38).push_back(72);
600 
601  fBranchingRatios.at(39).push_back(0.157258);
602  fBranchingRatios.at(39).push_back(0.403226);
603  fBranchingRatios.at(39).push_back(0.237903);
604  fBranchingRatios.at(39).push_back(0.201613);
605  fDecayTo.at(39).push_back(62);
606  fDecayTo.at(39).push_back(70);
607  fDecayTo.at(39).push_back(71);
608  fDecayTo.at(39).push_back(72);
609 
610  fBranchingRatios.at(40).push_back(0.0740741);
611  fBranchingRatios.at(40).push_back(0.925926);
612  fDecayTo.at(40).push_back(52);
613  fDecayTo.at(40).push_back(72);
614 
615  fBranchingRatios.at(41).push_back(0.0535714);
616  fBranchingRatios.at(41).push_back(0.35119);
617  fBranchingRatios.at(41).push_back(0.595238);
618  fDecayTo.at(41).push_back(67);
619  fDecayTo.at(41).push_back(68);
620  fDecayTo.at(41).push_back(70);
621 
622  fBranchingRatios.at(42).push_back(0.00816803);
623  fBranchingRatios.at(42).push_back(0.0583431);
624  fBranchingRatios.at(42).push_back(0.350058);
625  fBranchingRatios.at(42).push_back(0.583431);
626  fDecayTo.at(42).push_back(49);
627  fDecayTo.at(42).push_back(62);
628  fDecayTo.at(42).push_back(71);
629  fDecayTo.at(42).push_back(72);
630 
631  fBranchingRatios.at(43).push_back(0.0961538);
632  fBranchingRatios.at(43).push_back(0.423077);
633  fBranchingRatios.at(43).push_back(0.480769);
634  fDecayTo.at(43).push_back(66);
635  fDecayTo.at(43).push_back(67);
636  fDecayTo.at(43).push_back(68);
637 
638  fBranchingRatios.at(44).push_back(0.450549);
639  fBranchingRatios.at(44).push_back(0.549451);
640  fDecayTo.at(44).push_back(69);
641  fDecayTo.at(44).push_back(72);
642 
643  fBranchingRatios.at(45).push_back(0.469925);
644  fBranchingRatios.at(45).push_back(0.0836466);
645  fBranchingRatios.at(45).push_back(0.446429);
646  fDecayTo.at(45).push_back(61);
647  fDecayTo.at(45).push_back(65);
648  fDecayTo.at(45).push_back(72);
649 
650  fBranchingRatios.at(46).push_back(0.0408163);
651  fBranchingRatios.at(46).push_back(0.510204);
652  fBranchingRatios.at(46).push_back(0.44898);
653  fDecayTo.at(46).push_back(67);
654  fDecayTo.at(46).push_back(70);
655  fDecayTo.at(46).push_back(71);
656 
657  fBranchingRatios.at(47).push_back(1);
658  fDecayTo.at(47).push_back(72);
659 
660  fBranchingRatios.at(48).push_back(0.641026);
661  fBranchingRatios.at(48).push_back(0.358974);
662  fDecayTo.at(48).push_back(58);
663  fDecayTo.at(48).push_back(69);
664 
665  fBranchingRatios.at(49).push_back(0.0458015);
666  fBranchingRatios.at(49).push_back(0.954198);
667  fDecayTo.at(49).push_back(66);
668  fDecayTo.at(49).push_back(70);
669 
670  fBranchingRatios.at(50).push_back(0.401639);
671  fBranchingRatios.at(50).push_back(0.188525);
672  fBranchingRatios.at(50).push_back(0.409836);
673  fDecayTo.at(50).push_back(61);
674  fDecayTo.at(50).push_back(69);
675  fDecayTo.at(50).push_back(72);
676 
677  fBranchingRatios.at(51).push_back(0.0188679);
678  fBranchingRatios.at(51).push_back(0.173585);
679  fBranchingRatios.at(51).push_back(0.754717);
680  fBranchingRatios.at(51).push_back(0.0528302);
681  fDecayTo.at(51).push_back(61);
682  fDecayTo.at(51).push_back(67);
683  fDecayTo.at(51).push_back(71);
684  fDecayTo.at(51).push_back(72);
685 
686  fBranchingRatios.at(52).push_back(0.0100849);
687  fBranchingRatios.at(52).push_back(0.00796178);
688  fBranchingRatios.at(52).push_back(0.530786);
689  fBranchingRatios.at(52).push_back(0.451168);
690  fDecayTo.at(52).push_back(59);
691  fDecayTo.at(52).push_back(68);
692  fDecayTo.at(52).push_back(70);
693  fDecayTo.at(52).push_back(71);
694 
695  fBranchingRatios.at(53).push_back(0.0379902);
696  fBranchingRatios.at(53).push_back(0.0490196);
697  fBranchingRatios.at(53).push_back(0.612745);
698  fBranchingRatios.at(53).push_back(0.300245);
699  fDecayTo.at(53).push_back(67);
700  fDecayTo.at(53).push_back(70);
701  fDecayTo.at(53).push_back(71);
702  fDecayTo.at(53).push_back(72);
703 
704  fBranchingRatios.at(54).push_back(0.106557);
705  fBranchingRatios.at(54).push_back(0.819672);
706  fBranchingRatios.at(54).push_back(0.0737705);
707  fDecayTo.at(54).push_back(59);
708  fDecayTo.at(54).push_back(68);
709  fDecayTo.at(54).push_back(70);
710 
711  fBranchingRatios.at(55).push_back(0.699301);
712  fBranchingRatios.at(55).push_back(0.300699);
713  fDecayTo.at(55).push_back(64);
714  fDecayTo.at(55).push_back(70);
715 
716  fBranchingRatios.at(56).push_back(1);
717  fDecayTo.at(56).push_back(71);
718 
719  fBranchingRatios.at(57).push_back(1);
720  fDecayTo.at(57).push_back(72);
721 
722  fBranchingRatios.at(58).push_back(0.888099);
723  fBranchingRatios.at(58).push_back(0.111901);
724  fDecayTo.at(58).push_back(69);
725  fDecayTo.at(58).push_back(72);
726 
727  fBranchingRatios.at(59).push_back(0.00647298);
728  fBranchingRatios.at(59).push_back(0.752672);
729  fBranchingRatios.at(59).push_back(0.165588);
730  fBranchingRatios.at(59).push_back(0.0752672);
731  fDecayTo.at(59).push_back(65);
732  fDecayTo.at(59).push_back(70);
733  fDecayTo.at(59).push_back(71);
734  fDecayTo.at(59).push_back(72);
735 
736  fBranchingRatios.at(60).push_back(0.0708556);
737  fBranchingRatios.at(60).push_back(0.668449);
738  fBranchingRatios.at(60).push_back(0.260695);
739  fDecayTo.at(60).push_back(65);
740  fDecayTo.at(60).push_back(71);
741  fDecayTo.at(60).push_back(72);
742 
743  fBranchingRatios.at(61).push_back(0.166667);
744  fBranchingRatios.at(61).push_back(0.833333);
745  fDecayTo.at(61).push_back(69);
746  fDecayTo.at(61).push_back(72);
747 
748  fBranchingRatios.at(62).push_back(0.0898551);
749  fBranchingRatios.at(62).push_back(0.57971);
750  fBranchingRatios.at(62).push_back(0.330435);
751  fDecayTo.at(62).push_back(67);
752  fDecayTo.at(62).push_back(68);
753  fDecayTo.at(62).push_back(70);
754 
755  fBranchingRatios.at(63).push_back(0.813008);
756  fBranchingRatios.at(63).push_back(0.186992);
757  fDecayTo.at(63).push_back(71);
758  fDecayTo.at(63).push_back(72);
759 
760  fBranchingRatios.at(64).push_back(0.29078);
761  fBranchingRatios.at(64).push_back(0.70922);
762  fDecayTo.at(64).push_back(70);
763  fDecayTo.at(64).push_back(71);
764 
765  fBranchingRatios.at(65).push_back(0.05);
766  fBranchingRatios.at(65).push_back(0.08);
767  fBranchingRatios.at(65).push_back(0.5);
768  fBranchingRatios.at(65).push_back(0.37);
769  fDecayTo.at(65).push_back(69);
770  fDecayTo.at(65).push_back(70);
771  fDecayTo.at(65).push_back(71);
772  fDecayTo.at(65).push_back(72);
773 
774  fBranchingRatios.at(66).push_back(0.398406);
775  fBranchingRatios.at(66).push_back(0.310757);
776  fBranchingRatios.at(66).push_back(0.290837);
777  fDecayTo.at(66).push_back(70);
778  fDecayTo.at(66).push_back(71);
779  fDecayTo.at(66).push_back(72);
780 
781  fBranchingRatios.at(67).push_back(0.819672);
782  fBranchingRatios.at(67).push_back(0.180328);
783  fDecayTo.at(67).push_back(70);
784  fDecayTo.at(67).push_back(71);
785 
786  fBranchingRatios.at(68).push_back(0.186992);
787  fBranchingRatios.at(68).push_back(0.813008);
788  fDecayTo.at(68).push_back(70);
789  fDecayTo.at(68).push_back(71);
790 
791  fBranchingRatios.at(69).push_back(1);
792  fDecayTo.at(69).push_back(72);
793 
794  fBranchingRatios.at(70).push_back(1);
795  fDecayTo.at(70).push_back(71);
796 
797  fBranchingRatios.at(71).push_back(1);
798  fDecayTo.at(71).push_back(72);
799 
800  // Info from table VII (http://prc.aps.org/pdf/PRC/v58/i6/p3677_1)
801  // in MeV
802  double startEnergyLevels[] = { 2.281, 2.752, 2.937,
803  3.143, 3.334, 3.569, 3.652, 3.786, 3.861, 4.067, 4.111, 4.267,
804  4.364, 4.522, 4.655, 4.825, 5.017, 5.080, 5.223, 5.696, 6.006 };
805 
806  fStartEnergyLevels.insert(fStartEnergyLevels.end(), &startEnergyLevels[0],
807  &startEnergyLevels[fNumberOfStartLevels]);
808 
809  double b[] = { 0.9, 1.5, 0.11, 0.06, 0.04, 0.01,
810  0.16, 0.26, 0.01, 0.05, 0.11, 0.29,
811  3.84, 0.31, 0.38, 0.47, 0.36, 0.23,
812  0.03, 0.11, 0.13 };
813 
814  fB.insert(fB.end(), &b[0], &b[fNumberOfStartLevels]);
815 
816  double energyLevels[] = { 7.4724, 6.2270, 5.06347,
817  4.99294, 4.8756, 4.87255, 4.78865, 4.744093, 4.53706,
818  4.47299, 4.41936, 4.39588, 4.3840, 4.3656, 4.28052,
819  4.25362, 4.21307, 4.18003, 4.14901, 4.11084, 4.10446,
820  4.02035, 3.92390, 3.88792, 3.86866, 3.840228, 3.82143,
821  3.79757, 3.76779, 3.73848, 3.663739, 3.62995, 3.59924,
822  3.55697, 3.48621, 3.439144, 3.41434, 3.39363, 3.36803,
823  3.22867, 3.15381, 3.14644, 3.12836, 3.109721, 3.1002,
824  3.02795, 2.98587, 2.9508, 2.87901, 2.80788, 2.7874,
825  2.786644, 2.75672, 2.74691, 2.730372, 2.625990, 2.57593,
826  2.558, 2.54277, 2.419171, 2.397165, 2.290493, 2.289871,
827  2.26040, 2.103668, 2.069809, 2.047354, 1.959068, 1.643639,
828  0.891398, 0.8001427, 0.0298299, 0.0 };
829 
830  fEnergyLevels.insert(fEnergyLevels.end(), &energyLevels[0],
831  &energyLevels[fNumberOfLevels]);
832 
833  // I have to put checks that the arrays really
834  // have as many elements as they should
835 
836  }
837 
838  //----------------------------------------------------------------------------
839  // Simulate particles and fill truth
841  CLHEP::HepRandomEngine& engine) const
842  {
843 
844  bool success = false;
845  while (!success) {
846  double neutrinoEnergy = GetNeutrinoEnergy(engine);
847  double neutrinoTime = GetNeutrinoTime (engine);
848  success = ProcessOneNeutrino(truth, neutrinoEnergy, neutrinoTime, engine);
849  }
850 
851  }
852 
853  //----------------------------------------------------------------------------
854  // Simulate particles and fill truth
855  // Return true if one neutrino is processed, false otherwise
857  double neutrinoEnergy, double neutrinoTime,
858  CLHEP::HepRandomEngine& engine) const
859  {
860 
861  CLHEP::RandFlat randFlat(engine);
862 
863  int highestLevel = 0;
864  std::vector< double > levelCrossSections =
865  CalculateCrossSections(neutrinoEnergy, highestLevel);
866 
867  double totalCrossSection = 0;
868  // Calculating total cross section
869  for (std::vector< double >::iterator crossSection =
870  levelCrossSections.begin();
871  crossSection != levelCrossSections.end(); ++crossSection)
872  totalCrossSection += *crossSection;
873 
874  if (totalCrossSection == 0)
875  return false;
876 
877  std::vector< double > startLevelProbabilities;
878  // Calculating each level's probability
879  for (std::vector< double >::iterator crossSection =
880  levelCrossSections.begin();
881  crossSection != levelCrossSections.end(); ++crossSection)
882  startLevelProbabilities.push_back((*crossSection)/totalCrossSection);
883 
884  double randomNumber = randFlat.fire();
885  double tprob = 0;
886  int chosenStartLevel = -1;
887  // Picking a starting level
888  for (int level = 0; level < highestLevel; ++level)
889  {
890  if (randomNumber < (startLevelProbabilities.at(level) + tprob))
891  {
892  chosenStartLevel = level;
893  break;
894  }
895  tprob += startLevelProbabilities.at(level);
896  }
897 
898  int lastLevel = -1;
899  int level = -1;
900 
901  // Time delays
902  // Seems like this is not implemented yet
903  std::vector< double > levelDelay;
904  for (int n = 0; n < fNumberOfLevels; ++n) levelDelay.push_back(0.0);
905 
906  // The highest n for which fStartEnergyLevels.at(chosenStartLevel)
907  // is higher than fEnergyLevels.at(n)
908  int highestHigher = 0;
909  // The lowest n for which fStartEnergyLevels.at(chosenStartLevel)
910  // is lower than fEnergyLevels.at(n)
911  int lowestLower = 0;
912 
913  // Finding lowestLower and highestHigher
914  for (int n = 0; n < fNumberOfLevels; ++n)
915  {
916  if (fStartEnergyLevels.at(chosenStartLevel) < fEnergyLevels.at(n))
917  lowestLower = n;
918  if (fStartEnergyLevels.at(chosenStartLevel) > fEnergyLevels.at(n))
919  {
920  highestHigher = n;
921  break;
922  }
923  }
924 
925  if (std::abs(fStartEnergyLevels.at(chosenStartLevel) -
926  fEnergyLevels.at(lowestLower))
927  < std::abs(fStartEnergyLevels.at(chosenStartLevel)
928  - fEnergyLevels.at(highestHigher)))
929  {
930  lastLevel = lowestLower;
931  level = lowestLower;
932  }
933  // If the chosen start level energy is closest to the lowest level energy
934  // that it's lower than than the highest level energy
935  // that it's higher than, it starts at the level of
936  // the lowest level energy that it's lower than
937 
938  if (std::abs(fStartEnergyLevels.at(chosenStartLevel)
939  - fEnergyLevels.at(highestHigher))
940  < std::abs(fStartEnergyLevels.at(chosenStartLevel)
941  - fEnergyLevels.at(lowestLower)))
942  {
943  lastLevel = highestHigher;
944  level = highestHigher;
945  }
946  // If the chosen start level energy is closest to the highest level energy
947  // that it's higher than than the lowest level energy that it's lower than,
948  // it starts at the level of the highest level energy that it's higher than
949 
950  std::vector< double > vertex = GetUniformPosition(engine);
951 
952  std::vector< double > neutrinoDirection = GetIsotropicDirection(engine);
953 
954 
955 
956  //
957  // Add the first particle (the neutrino) to the truth list...
958  //
959 
960  // NOTE: all of the values below calculated for the neutrino should be checked!
961 
962  // Primary particles have negative IDs
963  int neutrinoTrackID = -1*(truth.NParticles() + 1);
964  std::string primary("primary");
965 
966  int nuePDG = 12;
967  double neutrinoP = neutrinoEnergy/1000.0; // use GeV...
968  double neutrinoPx = neutrinoDirection.at(0)*neutrinoP;
969  double neutrinoPy = neutrinoDirection.at(1)*neutrinoP;
970  double neutrinoPz = neutrinoDirection.at(2)*neutrinoP;
971 
972 
973 
974  // Create the neutrino:
975  // * set the mother to -1 since it is primary
976  // * set the mass to something < 0 so that the constructor looks up the mass from the pdg tables
977  // * set the status bit to 0 so that geant doesn't waste any CPU tracking the neutrino
978  simb::MCParticle neutrino(neutrinoTrackID, nuePDG, primary, -1, -1, 0);
979 
980  TLorentzVector neutrinoPosition(vertex.at(0), vertex.at(1),
981  vertex.at(2), neutrinoTime);
982  TLorentzVector neutrinoMomentum(neutrinoPx, neutrinoPy,
983  neutrinoPz, neutrinoEnergy/1000.0);
984  neutrino.AddTrajectoryPoint(neutrinoPosition, neutrinoMomentum);
985 
986  truth.Add(neutrino);
987 
988 
989 
990  //
991  // Adding the electron to truth
992  //
993 
994  // In MeV
995  double electronEnergy = neutrinoEnergy -
996  (fEnergyLevels.at(lastLevel) + 1.5);
997  double electronEnergyGeV = electronEnergy/1000.0;
998 
999  double electronM = 0.000511;
1000  double electronP = std::sqrt(std::pow(electronEnergyGeV,2)
1001  - std::pow(electronM,2));
1002 
1003  // For the moment, choose an isotropic direction for the electron.
1004  std::vector< double > electronDirection = GetIsotropicDirection(engine);
1005  double electronPx = electronDirection.at(0)*electronP;
1006  double electronPy = electronDirection.at(1)*electronP;
1007  double electronPz = electronDirection.at(2)*electronP;
1008 
1009  // Primary particles have negative IDs
1010  int trackID = -1*(truth.NParticles() + 1);
1011  int electronPDG = 11;
1012  simb::MCParticle electron(trackID, electronPDG, primary);
1013 
1014  TLorentzVector electronPosition(vertex.at(0), vertex.at(1),
1015  vertex.at(2), neutrinoTime);
1016  TLorentzVector electronMomentum(electronPx, electronPy,
1017  electronPz, electronEnergyGeV);
1018  electron.AddTrajectoryPoint(electronPosition, electronMomentum);
1019 
1020  truth.Add(electron);
1021 
1022  double ttime = neutrinoTime;
1023  int noMoreDecay = 0;
1024 
1025  // Level loop
1026  int groundLevel = fNumberOfLevels - 1;
1027  while (level != groundLevel)
1028  {
1029 
1030  double rl = randFlat.fire();
1031 
1032  int decayNum = 0;
1033 
1034  tprob = 0; // Used this variable above for cross section
1035 
1036  // Decay loop
1037  for (unsigned int iLevel = 0;
1038  iLevel < fBranchingRatios.at(level).size(); ++iLevel)
1039  {
1040 
1041  if (rl < (fBranchingRatios.at(level).at(iLevel) + tprob))
1042  {
1043 
1044  // We have a decay
1045 
1046  level = fDecayTo.at(level).at(decayNum);
1047 
1048  // Really should be using a map
1049  double gammaEnergy = fEnergyLevels.at(lastLevel) -
1050  fEnergyLevels.at(level);
1051  double gammaEnergyGeV = gammaEnergy/1000;
1052 
1053  std::vector< double > gammaDirection = GetIsotropicDirection(engine);
1054  double gammaPx = gammaDirection.at(0)*gammaEnergyGeV;
1055  double gammaPy = gammaDirection.at(1)*gammaEnergyGeV;
1056  double gammaPz = gammaDirection.at(2)*gammaEnergyGeV;
1057 
1058  //double gammaM = 0.0; // unused
1059 
1060  double gammaTime = (-TMath::Log(randFlat.fire())/
1061  (1/(levelDelay.at(lastLevel)))) + ttime;
1062 
1063  // Adding the gamma to truth
1064 
1065  // Primary particles have negative IDs
1066  trackID = -1*(truth.NParticles() + 1);
1067  int gammaPDG = 22;
1068  simb::MCParticle gamma(trackID, gammaPDG, primary); // NOTE: should the gammas be primary?
1069 
1070  TLorentzVector gammaPosition(vertex.at(0), vertex.at(1),
1071  vertex.at(2), neutrinoTime);
1072  TLorentzVector gammaMomentum(gammaPx, gammaPy,
1073  gammaPz, gammaEnergyGeV);
1074  gamma.AddTrajectoryPoint(gammaPosition, gammaMomentum);
1075 
1076  truth.Add(gamma);
1077 
1078  lastLevel = level;
1079  ttime = gammaTime;
1080 
1081  break;
1082 
1083  }
1084 
1085  if ((tprob + fBranchingRatios.at(level).at(iLevel)) > 1) {
1086  std::cout << "(tprob + *iLevel) > 1" << std::endl;
1087  noMoreDecay = 1; // If it doesn't do any more gamma decay
1088  break;
1089  }
1090 
1091  ++decayNum;
1092  tprob += fBranchingRatios.at(level).at(iLevel);
1093 
1094  } // End of decay loop
1095 
1096  if (noMoreDecay == 1) break;
1097 
1098  } // End of level loop
1099 
1100  if (level != 72)
1101  {
1102  std::cout << "level != 72" << std::endl;
1103  std::cout << "level = " << level << std::endl;
1104  }
1105 
1106  // Set the neutrino for the MCTruth object:
1107  // NOTE: currently these parameters are all pretty much a guess...
1108  truth.SetNeutrino(simb::kCC,
1109  simb::kQE, // not sure what the mode should be here, assumed that these are all QE???
1110  simb::kCCQE, // not sure what the int_type should be either...
1111  0, // target is AR40? not sure how to specify that...
1112  0, // nucleon pdg ???
1113  0, // quark pdg ???
1114  -1.0, // W ??? - not sure how to calculate this from the above
1115  -1.0, // X ??? - not sure how to calculate this from the above
1116  -1.0, // Y ??? - not sure how to calculate this from the above
1117  -1.0); // Qsqr ??? - not sure how to calculate this from the above
1118 
1119  return true;
1120 
1121  }
1122 
1123  //----------------------------------------------------------------------------
1124  // Calculate cross sections for neutrino with neutrinoEnergy to excite
1125  // the final nucleus with highestLevel being the highest possible level.
1126  // highestLevel is output, so, whatever integer was used as the argument,
1127  // it may have a different value after this function is executed
1128  std::vector< double > NueAr40CCGenerator::CalculateCrossSections
1129  (double neutrinoEnergy, int& highestLevel) const
1130  {
1131 
1132  highestLevel = 0;
1133  std::vector< double > levelCrossSections;
1134 
1135  // Loop through energy levels, if neutrino has enough energy,
1136  // calculate cross section
1137  for (int level = 0; level < fNumberOfStartLevels; ++level)
1138  {
1139  // Electron energy in keV
1140  double w = (neutrinoEnergy - (fStartEnergyLevels.at(level) + 1.5))*1000;
1141 
1142  double sigma = 0.0;
1143  if (neutrinoEnergy > (fStartEnergyLevels.at(level) + 1.5) && w >= 511.)
1144  {
1145  ++highestLevel;
1146  for (int n = 0; n <= level; ++n)
1147  {
1148  // Electron energy for level n to use in cross section
1149  w = (neutrinoEnergy - (fStartEnergyLevels.at(n) + 1.5))*1000;
1150  // Electron momentum in keV/c
1151  double p = std::sqrt(pow(w, 2) - pow(511.0, 2));
1152  // Fermi function approximation
1153  double f = std::sqrt(3.0634 + (0.6814/(w - 1)));
1154  // In cm^2*10^-42
1155  sigma += 1.6e-8*(p*w*f*fB.at(n));
1156  }
1157  }
1158  levelCrossSections.push_back(sigma);
1159  }
1160 
1161  return levelCrossSections;
1162 
1163  }
1164 
1165 }
intermediate_table::iterator iterator
std::vector< double > fStartEnergyLevels
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:257
bool ProcessOneNeutrino(simb::MCTruth &truth, double neutrinoEnergy, double neutrinoTime, CLHEP::HepRandomEngine &engine) const
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
std::string string
Definition: nybbler.cc:12
std::vector< double > GetUniformPosition(CLHEP::HepRandomEngine &engine) const
std::vector< double > fEnergyLevels
constexpr T pow(T x)
Definition: pow.h:72
STL namespace.
intermediate_table::const_iterator const_iterator
void CreateKinematicsVector(simb::MCTruth &truth, CLHEP::HepRandomEngine &engine) const
int NParticles() const
Definition: MCTruth.h:75
void SetNeutrino(int CCNC, int mode, int interactionType, int target, int nucleon, int quark, double w, double x, double y, double qsqr)
Definition: MCTruth.cxx:30
charged current quasi-elastic
Definition: MCNeutrino.h:96
std::vector< std::vector< double > > fBranchingRatios
Particle class.
std::vector< double > GetIsotropicDirection(CLHEP::HepRandomEngine &engine) const
int const nuePDG
Definition: includeROOT.h:153
int const gammaPDG
Definition: includeROOT.h:165
T abs(T value)
double GetNeutrinoTime(CLHEP::HepRandomEngine &engine) const
std::vector< double > CalculateCrossSections(double neutrinoEnergy, int &highestLevel) const
int GetNumberOfNeutrinos(CLHEP::HepRandomEngine &engine) const
std::void_t< T > n
NueAr40CCGenerator(fhicl::ParameterSet const &parameterSet)
T get(std::string const &key) const
Definition: ParameterSet.h:271
p
Definition: test.py:223
std::vector< double > fB
std::map< double, double > fEnergyProbabilityMap
double gamma(double KE, const simb::MCParticle *part)
double GetNeutrinoEnergy(CLHEP::HepRandomEngine &engine) const
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
Supernova neutrinos.
Definition: MCTruth.h:25
std::string find_file(std::string const &filename) const
Definition: search_path.cc:96
static bool * b
Definition: config.cpp:1043
std::vector< simb::MCTruth > Generate(CLHEP::HepRandomEngine &engine)
Event generator information.
Definition: MCTruth.h:32
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
std::vector< std::vector< double > > fActiveVolume
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
int const electronPDG
Definition: includeROOT.h:151
std::vector< std::vector< int > > fDecayTo
Event Generation using GENIE, cosmics or single particles.
int bool
Definition: qglobal.h:345
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)
vertex reconstruction