Public Member Functions | Private Member Functions | Private Attributes | List of all members
evgen::NueAr40CCGenerator Class Reference

#include <NueAr40CCGenerator.h>

Public Member Functions

 NueAr40CCGenerator (fhicl::ParameterSet const &parameterSet)
 
std::vector< simb::MCTruth > Generate (CLHEP::HepRandomEngine &engine)
 

Private Member Functions

std::vector< double > GetIsotropicDirection (CLHEP::HepRandomEngine &engine) const
 
std::vector< double > GetUniformPosition (CLHEP::HepRandomEngine &engine) const
 
int GetNumberOfNeutrinos (CLHEP::HepRandomEngine &engine) const
 
double GetNeutrinoTime (CLHEP::HepRandomEngine &engine) const
 
double GetNeutrinoEnergy (CLHEP::HepRandomEngine &engine) const
 
void ReadNeutrinoSpectrum ()
 
void InitializeVectors ()
 
void CreateKinematicsVector (simb::MCTruth &truth, CLHEP::HepRandomEngine &engine) const
 
bool ProcessOneNeutrino (simb::MCTruth &truth, double neutrinoEnergy, double neutrinoTime, CLHEP::HepRandomEngine &engine) const
 
std::vector< double > CalculateCrossSections (double neutrinoEnergy, int &highestLevel) const
 

Private Attributes

std::map< double, double > fEnergyProbabilityMap
 
int fNumberOfLevels
 
int fNumberOfStartLevels
 
std::vector< std::vector< double > > fBranchingRatios
 
std::vector< std::vector< int > > fDecayTo
 
std::vector< double > fStartEnergyLevels
 
std::vector< double > fB
 
std::vector< double > fEnergyLevels
 
bool fMonoenergeticNeutrinos
 
double fNeutrinoEnergy
 
std::string fEnergySpectrumFileName
 
bool fUsePoissonDistribution
 
bool fAllowZeroNeutrinos
 
int fNumberOfNeutrinos
 
double fNeutrinoTimeBegin
 
double fNeutrinoTimeEnd
 
std::vector< std::vector< double > > fActiveVolume
 

Detailed Description

Definition at line 26 of file NueAr40CCGenerator.h.

Constructor & Destructor Documentation

evgen::NueAr40CCGenerator::NueAr40CCGenerator ( fhicl::ParameterSet const &  parameterSet)

Definition at line 35 of file NueAr40CCGenerator.cxx.

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 
64  if (!fMonoenergeticNeutrinos) ReadNeutrinoSpectrum();
65 
66  InitializeVectors();
67 
68  }
string
std::string string
Definition: nybbler.cc:12
evgen::NueAr40CCGenerator::fBranchingRatios
std::vector< std::vector< double > > fBranchingRatios
Definition: NueAr40CCGenerator.h:83
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

Member Function Documentation

std::vector< double > evgen::NueAr40CCGenerator::CalculateCrossSections ( double  neutrinoEnergy,
int &  highestLevel 
) const
private

Definition at line 1129 of file NueAr40CCGenerator.cxx.

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  }
evgen::NueAr40CCGenerator::fStartEnergyLevels
std::vector< double > fStartEnergyLevels
Definition: NueAr40CCGenerator.h:85
keras_to_tensorflow.f
f
Definition: keras_to_tensorflow.py:162
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
cet::n
std::void_t< T > n
Definition: metaprogramming.h:28
generate_datataset.level
level
Definition: generate_datataset.py:30
test.p
p
Definition: test.py:223
evgen::NueAr40CCGenerator::fB
std::vector< double > fB
Definition: NueAr40CCGenerator.h:86
void evgen::NueAr40CCGenerator::CreateKinematicsVector ( simb::MCTruth &  truth,
CLHEP::HepRandomEngine &  engine 
) const
private

Definition at line 840 of file NueAr40CCGenerator.cxx.

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  }
evgen::NueAr40CCGenerator::ProcessOneNeutrino
bool ProcessOneNeutrino(simb::MCTruth &truth, double neutrinoEnergy, double neutrinoTime, CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:856
evgen::NueAr40CCGenerator::GetNeutrinoTime
double GetNeutrinoTime(CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:161
evgen::NueAr40CCGenerator::GetNeutrinoEnergy
double GetNeutrinoEnergy(CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:174
std::vector< simb::MCTruth > evgen::NueAr40CCGenerator::Generate ( CLHEP::HepRandomEngine &  engine)

Definition at line 72 of file NueAr40CCGenerator.cxx.

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;
81  truth.SetOrigin(simb::kSuperNovaNeutrino);
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  }
simb::MCTruth::SetOrigin
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
evgen::NueAr40CCGenerator::CreateKinematicsVector
void CreateKinematicsVector(simb::MCTruth &truth, CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:840
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:75
evgen::NueAr40CCGenerator::GetNumberOfNeutrinos
int GetNumberOfNeutrinos(CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:143
simb::kSuperNovaNeutrino
Supernova neutrinos.
Definition: MCTruth.h:25
simb::MCTruth
Event generator information.
Definition: MCTruth.h:32
std::vector< double > evgen::NueAr40CCGenerator::GetIsotropicDirection ( CLHEP::HepRandomEngine &  engine) const
private

Definition at line 99 of file NueAr40CCGenerator.cxx.

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  }
truth_ana.theta
theta
Definition: truth_ana.py:138
double evgen::NueAr40CCGenerator::GetNeutrinoEnergy ( CLHEP::HepRandomEngine &  engine) const
private

Definition at line 174 of file NueAr40CCGenerator.cxx.

175  {
176 
177  if (fMonoenergeticNeutrinos) return fNeutrinoEnergy;
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  }
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
evgen::NueAr40CCGenerator::fEnergyProbabilityMap
std::map< double, double > fEnergyProbabilityMap
Definition: NueAr40CCGenerator.h:77
util::quantities::second
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
double evgen::NueAr40CCGenerator::GetNeutrinoTime ( CLHEP::HepRandomEngine &  engine) const
private

Definition at line 161 of file NueAr40CCGenerator.cxx.

162  {
163 
164  CLHEP::RandFlat randFlat(engine);
165 
166  return randFlat.fire(fNeutrinoTimeBegin, fNeutrinoTimeEnd);
167 
168  }
int evgen::NueAr40CCGenerator::GetNumberOfNeutrinos ( CLHEP::HepRandomEngine &  engine) const
private

Definition at line 143 of file NueAr40CCGenerator.cxx.

144  {
145 
146  if (fUsePoissonDistribution)
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  }
std::vector< double > evgen::NueAr40CCGenerator::GetUniformPosition ( CLHEP::HepRandomEngine &  engine) const
private

Definition at line 122 of file NueAr40CCGenerator.cxx.

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  }
evgen::NueAr40CCGenerator::fActiveVolume
std::vector< std::vector< double > > fActiveVolume
Definition: NueAr40CCGenerator.h:114
void evgen::NueAr40CCGenerator::InitializeVectors ( )
private

Definition at line 258 of file NueAr40CCGenerator.cxx.

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  }
evgen::NueAr40CCGenerator::fStartEnergyLevels
std::vector< double > fStartEnergyLevels
Definition: NueAr40CCGenerator.h:85
evgen::NueAr40CCGenerator::fEnergyLevels
std::vector< double > fEnergyLevels
Definition: NueAr40CCGenerator.h:87
evgen::NueAr40CCGenerator::fBranchingRatios
std::vector< std::vector< double > > fBranchingRatios
Definition: NueAr40CCGenerator.h:83
evgen::NueAr40CCGenerator::fB
std::vector< double > fB
Definition: NueAr40CCGenerator.h:86
b
static bool * b
Definition: config.cpp:1043
evgen::NueAr40CCGenerator::fDecayTo
std::vector< std::vector< int > > fDecayTo
Definition: NueAr40CCGenerator.h:84
bool evgen::NueAr40CCGenerator::ProcessOneNeutrino ( simb::MCTruth &  truth,
double  neutrinoEnergy,
double  neutrinoTime,
CLHEP::HepRandomEngine &  engine 
) const
private

Definition at line 856 of file NueAr40CCGenerator.cxx.

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  }
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
evgen::NueAr40CCGenerator::fStartEnergyLevels
std::vector< double > fStartEnergyLevels
Definition: NueAr40CCGenerator.h:85
string
std::string string
Definition: nybbler.cc:12
evgen::NueAr40CCGenerator::GetUniformPosition
std::vector< double > GetUniformPosition(CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:122
evgen::NueAr40CCGenerator::fEnergyLevels
std::vector< double > fEnergyLevels
Definition: NueAr40CCGenerator.h:87
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:75
simb::MCTruth::SetNeutrino
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
simb::kCCQE
charged current quasi-elastic
Definition: MCNeutrino.h:96
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:99
nuePDG
int const nuePDG
Definition: includeROOT.h:153
gammaPDG
int const gammaPDG
Definition: includeROOT.h:165
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
evgen::NueAr40CCGenerator::CalculateCrossSections
std::vector< double > CalculateCrossSections(double neutrinoEnergy, int &highestLevel) const
Definition: NueAr40CCGenerator.cxx:1129
cet::n
std::void_t< T > n
Definition: metaprogramming.h:28
generate_datataset.level
level
Definition: generate_datataset.py:30
pionana::gamma
double gamma(double KE, const simb::MCParticle *part)
Definition: TruthAnalyzer_module.cc:55
simb::MCTruth::Add
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
electronPDG
int const electronPDG
Definition: includeROOT.h:151
evgen::NueAr40CCGenerator::fDecayTo
std::vector< std::vector< int > > fDecayTo
Definition: NueAr40CCGenerator.h:84
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
vertex
vertex reconstruction
Definition: AggregateVertexAna_module.cc:32
void evgen::NueAr40CCGenerator::ReadNeutrinoSpectrum ( )
private

Definition at line 217 of file NueAr40CCGenerator.cxx.

218  {
219 
220  cet::search_path searchPath("FW_SEARCH_PATH");
221  std::string directoryName = "SupernovaNeutrinos/" +
222  fEnergySpectrumFileName;
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  }
string
std::string string
Definition: nybbler.cc:12
evgen::NueAr40CCGenerator::fEnergyProbabilityMap
std::map< double, double > fEnergyProbabilityMap
Definition: NueAr40CCGenerator.h:77
cet::search_path
Definition: search_path.h:36
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Data Documentation

std::vector< std::vector < double > > evgen::NueAr40CCGenerator::fActiveVolume
private

Definition at line 114 of file NueAr40CCGenerator.h.

bool evgen::NueAr40CCGenerator::fAllowZeroNeutrinos
private

Definition at line 103 of file NueAr40CCGenerator.h.

std::vector< double > evgen::NueAr40CCGenerator::fB
private

Definition at line 86 of file NueAr40CCGenerator.h.

std::vector< std::vector< double > > evgen::NueAr40CCGenerator::fBranchingRatios
private

Definition at line 83 of file NueAr40CCGenerator.h.

std::vector< std::vector< int > > evgen::NueAr40CCGenerator::fDecayTo
private

Definition at line 84 of file NueAr40CCGenerator.h.

std::vector< double > evgen::NueAr40CCGenerator::fEnergyLevels
private

Definition at line 87 of file NueAr40CCGenerator.h.

std::map< double, double > evgen::NueAr40CCGenerator::fEnergyProbabilityMap
private

Definition at line 77 of file NueAr40CCGenerator.h.

std::string evgen::NueAr40CCGenerator::fEnergySpectrumFileName
private

Definition at line 95 of file NueAr40CCGenerator.h.

bool evgen::NueAr40CCGenerator::fMonoenergeticNeutrinos
private

Definition at line 90 of file NueAr40CCGenerator.h.

double evgen::NueAr40CCGenerator::fNeutrinoEnergy
private

Definition at line 92 of file NueAr40CCGenerator.h.

double evgen::NueAr40CCGenerator::fNeutrinoTimeBegin
private

Definition at line 109 of file NueAr40CCGenerator.h.

double evgen::NueAr40CCGenerator::fNeutrinoTimeEnd
private

Definition at line 110 of file NueAr40CCGenerator.h.

int evgen::NueAr40CCGenerator::fNumberOfLevels
private

Definition at line 79 of file NueAr40CCGenerator.h.

int evgen::NueAr40CCGenerator::fNumberOfNeutrinos
private

Definition at line 106 of file NueAr40CCGenerator.h.

int evgen::NueAr40CCGenerator::fNumberOfStartLevels
private

Definition at line 80 of file NueAr40CCGenerator.h.

std::vector< double > evgen::NueAr40CCGenerator::fStartEnergyLevels
private

Definition at line 85 of file NueAr40CCGenerator.h.

bool evgen::NueAr40CCGenerator::fUsePoissonDistribution
private

Definition at line 99 of file NueAr40CCGenerator.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11