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