LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
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 27 of file NueAr40CCGenerator.h.

Constructor & Destructor Documentation

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

Definition at line 33 of file NueAr40CCGenerator.cxx.

References fActiveVolume, fMonoenergeticNeutrinos, fhicl::ParameterSet::get(), InitializeVectors(), and ReadNeutrinoSpectrum().

35  : fNumberOfLevels ( 73 )
36  , fNumberOfStartLevels ( 21 )
37  , fBranchingRatios (fNumberOfLevels )
38  , fDecayTo (fNumberOfLevels )
39  , fMonoenergeticNeutrinos
40  (parameterSet.get< bool >("MonoenergeticNeutrinos"))
41  , fNeutrinoEnergy
42  (parameterSet.get< double >("NeutrinoEnergy" ))
43  , fEnergySpectrumFileName
44  (parameterSet.get< std::string >("EnergySpectrumFileName"))
45  , fUsePoissonDistribution
46  (parameterSet.get< bool >("UsePoissonDistribution"))
47  , fAllowZeroNeutrinos
48  (parameterSet.get< bool >("AllowZeroNeutrinos" ))
49  , fNumberOfNeutrinos
50  (parameterSet.get< int >("NumberOfNeutrinos" ))
51  , fNeutrinoTimeBegin
52  (parameterSet.get< double >("NeutrinoTimeBegin" ))
53  , fNeutrinoTimeEnd
54  (parameterSet.get< double >("NeutrinoTimeEnd" ))
55  {
56 
57  fActiveVolume.push_back
58  (parameterSet.get< std::vector< double > >("ActiveVolume0"));
59  fActiveVolume.push_back
60  (parameterSet.get< std::vector< double > >("ActiveVolume1"));
61 
62  if (!fMonoenergeticNeutrinos) ReadNeutrinoSpectrum();
63 
64  InitializeVectors();
65 
66  }
evgen::NueAr40CCGenerator::fBranchingRatios
std::vector< std::vector< double > > fBranchingRatios
Definition: NueAr40CCGenerator.h:84
evgen::NueAr40CCGenerator::fActiveVolume
std::vector< std::vector< double > > fActiveVolume
Definition: NueAr40CCGenerator.h:115
evgen::NueAr40CCGenerator::fDecayTo
std::vector< std::vector< int > > fDecayTo
Definition: NueAr40CCGenerator.h:85

Member Function Documentation

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

Definition at line 1127 of file NueAr40CCGenerator.cxx.

References f, fB, fNumberOfStartLevels, fStartEnergyLevels, n, and w.

Referenced by ProcessOneNeutrino().

1128  {
1129 
1130  highestLevel = 0;
1131  std::vector< double > levelCrossSections;
1132 
1133  // Loop through energy levels, if neutrino has enough energy,
1134  // calculate cross section
1135  for (int level = 0; level < fNumberOfStartLevels; ++level)
1136  {
1137  // Electron energy in keV
1138  double w = (neutrinoEnergy - (fStartEnergyLevels.at(level) + 1.5))*1000;
1139 
1140  double sigma = 0.0;
1141  if (neutrinoEnergy > (fStartEnergyLevels.at(level) + 1.5) && w >= 511.)
1142  {
1143  ++highestLevel;
1144  for (int n = 0; n <= level; ++n)
1145  {
1146  // Electron energy for level n to use in cross section
1147  w = (neutrinoEnergy - (fStartEnergyLevels.at(n) + 1.5))*1000;
1148  // Electron momentum in keV/c
1149  double p = std::sqrt(pow(w, 2) - pow(511.0, 2));
1150  // Fermi function approximation
1151  double f = std::sqrt(3.0634 + (0.6814/(w - 1)));
1152  // In cm^2*10^-42
1153  sigma += 1.6e-8*(p*w*f*fB.at(n));
1154  }
1155  }
1156  levelCrossSections.push_back(sigma);
1157  }
1158 
1159  return levelCrossSections;
1160 
1161  }
evgen::NueAr40CCGenerator::fStartEnergyLevels
std::vector< double > fStartEnergyLevels
Definition: NueAr40CCGenerator.h:86
f
TFile f
Definition: plotHisto.C:6
evgen::NueAr40CCGenerator::fB
std::vector< double > fB
Definition: NueAr40CCGenerator.h:87
n
Char_t n[5]
Definition: comparison_ascii.C:55
w
Float_t w
Definition: plot.C:23
void evgen::NueAr40CCGenerator::CreateKinematicsVector ( simb::MCTruth &  truth,
CLHEP::HepRandomEngine &  engine 
) const
private

Definition at line 838 of file NueAr40CCGenerator.cxx.

References GetNeutrinoEnergy(), GetNeutrinoTime(), and ProcessOneNeutrino().

Referenced by Generate().

840  {
841 
842  bool success = false;
843  while (!success) {
844  double neutrinoEnergy = GetNeutrinoEnergy(engine);
845  double neutrinoTime = GetNeutrinoTime (engine);
846  success = ProcessOneNeutrino(truth, neutrinoEnergy, neutrinoTime, engine);
847  }
848 
849  }
evgen::NueAr40CCGenerator::ProcessOneNeutrino
bool ProcessOneNeutrino(simb::MCTruth &truth, double neutrinoEnergy, double neutrinoTime, CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:854
evgen::NueAr40CCGenerator::GetNeutrinoTime
double GetNeutrinoTime(CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:159
evgen::NueAr40CCGenerator::GetNeutrinoEnergy
double GetNeutrinoEnergy(CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:172
std::vector< simb::MCTruth > evgen::NueAr40CCGenerator::Generate ( CLHEP::HepRandomEngine &  engine)

Definition at line 70 of file NueAr40CCGenerator.cxx.

References CreateKinematicsVector(), GetIsotropicDirection(), GetNumberOfNeutrinos(), simb::kSuperNovaNeutrino, simb::MCTruth::NParticles(), and simb::MCTruth::SetOrigin().

Referenced by evgen::SNNueAr40CCGen::produce().

71  {
72 
73  std::vector<simb::MCTruth> truths;
74 
75  int NumberOfNu = this->GetNumberOfNeutrinos(engine);
76  for(int i = 0; i < NumberOfNu; ++i) {
77 
78  simb::MCTruth truth;
79  truth.SetOrigin(simb::kSuperNovaNeutrino);
80 
81  // Loop until at least one neutrino is simulated
82  while (!truth.NParticles()) {
83  CreateKinematicsVector(truth, engine);
84  }
85 
86  truths.push_back(truth);
87 
88  }
89 
90  return truths;
91 
92  }
simb::MCTruth::SetOrigin
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:78
evgen::NueAr40CCGenerator::CreateKinematicsVector
void CreateKinematicsVector(simb::MCTruth &truth, CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:838
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:72
evgen::NueAr40CCGenerator::GetNumberOfNeutrinos
int GetNumberOfNeutrinos(CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:141
simb::kSuperNovaNeutrino
Supernova neutrinos.
Definition: MCTruth.h:23
simb::MCTruth
Event generator information.
Definition: MCTruth.h:30
std::vector< double > evgen::NueAr40CCGenerator::GetIsotropicDirection ( CLHEP::HepRandomEngine &  engine) const
private

Definition at line 97 of file NueAr40CCGenerator.cxx.

References GetUniformPosition().

Referenced by Generate(), and ProcessOneNeutrino().

98  {
99 
100  CLHEP::RandFlat randFlat(engine);
101 
102  std::vector< double > isotropicDirection;
103 
104  double phi = 2*TMath::Pi()*randFlat.fire();
105  double cosTheta = 2*randFlat.fire() - 1;
106  double theta = TMath::ACos(cosTheta);
107 
108  // x, y, z
109  isotropicDirection.push_back(cos(phi)*sin(theta));
110  isotropicDirection.push_back(sin(phi)*sin(theta));
111  isotropicDirection.push_back(cosTheta);
112 
113  return isotropicDirection;
114 
115  }
double evgen::NueAr40CCGenerator::GetNeutrinoEnergy ( CLHEP::HepRandomEngine &  engine) const
private

Definition at line 172 of file NueAr40CCGenerator.cxx.

References fEnergyProbabilityMap, fMonoenergeticNeutrinos, and fNeutrinoEnergy.

Referenced by CreateKinematicsVector(), and GetNeutrinoTime().

173  {
174 
175  if (fMonoenergeticNeutrinos) return fNeutrinoEnergy;
176 
177  CLHEP::RandFlat randFlat(engine);
178 
179  double neutrinoEnergy = 0.0;
180 
181  double randomNumber = randFlat.fire();
182 
183  // We need this to get a previous entry in the map
184  std::pair< double, double > previousPair;
185 
186  for (std::map< double, double >::const_iterator energyProbability =
187  fEnergyProbabilityMap.begin(); energyProbability !=
188  fEnergyProbabilityMap.end(); ++energyProbability)
189  {
190  if (randomNumber < energyProbability->second)
191  {
192  if (energyProbability != fEnergyProbabilityMap.begin())
193  {
194  neutrinoEnergy = energyProbability->first -
195  (energyProbability->second - randomNumber)*
196  (energyProbability->first - previousPair.first)/
197  (energyProbability->second - previousPair.second);
198  break;
199  }
200  else
201  {
202  neutrinoEnergy = energyProbability->first;
203  break;
204  }
205  }
206  previousPair = *energyProbability;
207  }
208 
209  return neutrinoEnergy;
210 
211  }
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:25
evgen::NueAr40CCGenerator::fEnergyProbabilityMap
std::map< double, double > fEnergyProbabilityMap
Definition: NueAr40CCGenerator.h:78
double evgen::NueAr40CCGenerator::GetNeutrinoTime ( CLHEP::HepRandomEngine &  engine) const
private

Definition at line 159 of file NueAr40CCGenerator.cxx.

References fNeutrinoTimeBegin, fNeutrinoTimeEnd, and GetNeutrinoEnergy().

Referenced by CreateKinematicsVector(), and GetNumberOfNeutrinos().

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

Definition at line 141 of file NueAr40CCGenerator.cxx.

References fAllowZeroNeutrinos, fNumberOfNeutrinos, fUsePoissonDistribution, and GetNeutrinoTime().

Referenced by Generate(), and GetUniformPosition().

142  {
143 
144  if (fUsePoissonDistribution)
145  {
146  CLHEP::RandPoisson randPoisson(engine);
147  int N = randPoisson.fire(fNumberOfNeutrinos);
148  if(N == 0 && !fAllowZeroNeutrinos) N = 1;
149  return N;
150  }
151 
152  return fNumberOfNeutrinos;
153 
154  }
std::vector< double > evgen::NueAr40CCGenerator::GetUniformPosition ( CLHEP::HepRandomEngine &  engine) const
private

Definition at line 120 of file NueAr40CCGenerator.cxx.

References fActiveVolume, and GetNumberOfNeutrinos().

Referenced by GetIsotropicDirection(), and ProcessOneNeutrino().

121  {
122 
123  CLHEP::RandFlat randFlat(engine);
124 
125  std::vector< double > position;
126 
127  position.push_back(randFlat.
128  fire(fActiveVolume.at(0).at(0), fActiveVolume.at(1).at(0)));
129  position.push_back(randFlat.
130  fire(fActiveVolume.at(0).at(1), fActiveVolume.at(1).at(1)));
131  position.push_back(randFlat.
132  fire(fActiveVolume.at(0).at(2), fActiveVolume.at(1).at(2)));
133 
134  return position;
135 
136  }
evgen::NueAr40CCGenerator::fActiveVolume
std::vector< std::vector< double > > fActiveVolume
Definition: NueAr40CCGenerator.h:115
void evgen::NueAr40CCGenerator::InitializeVectors ( )
private

Definition at line 256 of file NueAr40CCGenerator.cxx.

References fB, fBranchingRatios, fDecayTo, fEnergyLevels, fNumberOfLevels, fNumberOfStartLevels, and fStartEnergyLevels.

Referenced by NueAr40CCGenerator().

257  {
258 
259  // The level data -- there's probably a better way to initialize
260  fBranchingRatios.at(0).push_back(1.00);
261  fDecayTo.at(0).push_back(1);
262 
263  fBranchingRatios.at(1).push_back(1.00);
264  fDecayTo.at(1).push_back(4);
265 
266  fBranchingRatios.at(2).push_back(0.133891);
267  fBranchingRatios.at(2).push_back(0.41841);
268  fBranchingRatios.at(2).push_back(0.351464);
269  fBranchingRatios.at(2).push_back(0.0962343);
270  fDecayTo.at(2).push_back(38);
271  fDecayTo.at(2).push_back(43);
272  fDecayTo.at(2).push_back(59);
273  fDecayTo.at(2).push_back(72);
274 
275  fBranchingRatios.at(3).push_back(0.391705);
276  fBranchingRatios.at(3).push_back(0.460829);
277  fBranchingRatios.at(3).push_back(0.147465);
278  fDecayTo.at(3).push_back(61);
279  fDecayTo.at(3).push_back(65);
280  fDecayTo.at(3).push_back(71);
281 
282  fBranchingRatios.at(4).push_back(0.358974);
283  fBranchingRatios.at(4).push_back(0.641026);
284  fDecayTo.at(4).push_back(13);
285  fDecayTo.at(4).push_back(58);
286 
287  fBranchingRatios.at(5).push_back(0.247808);
288  fBranchingRatios.at(5).push_back(0.0468929);
289  fBranchingRatios.at(5).push_back(0.213496);
290  fBranchingRatios.at(5).push_back(0.11056);
291  fBranchingRatios.at(5).push_back(0.381243);
292  fDecayTo.at(5).push_back(40);
293  fDecayTo.at(5).push_back(52);
294  fDecayTo.at(5).push_back(67);
295  fDecayTo.at(5).push_back(71);
296  fDecayTo.at(5).push_back(72);
297 
298  fBranchingRatios.at(6).push_back(0.361025);
299  fBranchingRatios.at(6).push_back(0.056677);
300  fBranchingRatios.at(6).push_back(0.194099);
301  fBranchingRatios.at(6).push_back(0.388199);
302  fDecayTo.at(6).push_back(52);
303  fDecayTo.at(6).push_back(55);
304  fDecayTo.at(6).push_back(63);
305  fDecayTo.at(6).push_back(68);
306 
307  fBranchingRatios.at(7).push_back(0.0300963);
308  fBranchingRatios.at(7).push_back(0.0613965);
309  fBranchingRatios.at(7).push_back(0.0152488);
310  fBranchingRatios.at(7).push_back(0.0321027);
311  fBranchingRatios.at(7).push_back(0.0513644);
312  fBranchingRatios.at(7).push_back(0.0682183);
313  fBranchingRatios.at(7).push_back(0.073435);
314  fBranchingRatios.at(7).push_back(0.0100321);
315  fBranchingRatios.at(7).push_back(0.0120385);
316  fBranchingRatios.at(7).push_back(0.0922953);
317  fBranchingRatios.at(7).push_back(0.152488);
318  fBranchingRatios.at(7).push_back(0.401284);
319  fDecayTo.at(7).push_back(17);
320  fDecayTo.at(7).push_back(25);
321  fDecayTo.at(7).push_back(27);
322  fDecayTo.at(7).push_back(32);
323  fDecayTo.at(7).push_back(49);
324  fDecayTo.at(7).push_back(54);
325  fDecayTo.at(7).push_back(56);
326  fDecayTo.at(7).push_back(62);
327  fDecayTo.at(7).push_back(63);
328  fDecayTo.at(7).push_back(67);
329  fDecayTo.at(7).push_back(68);
330  fDecayTo.at(7).push_back(70);
331 
332  fBranchingRatios.at(8).push_back(0.0089877);
333  fBranchingRatios.at(8).push_back(0.06386);
334  fBranchingRatios.at(8).push_back(0.13245);
335  fBranchingRatios.at(8).push_back(0.473037);
336  fBranchingRatios.at(8).push_back(0.321665);
337  fDecayTo.at(8).push_back(43);
338  fDecayTo.at(8).push_back(56);
339  fDecayTo.at(8).push_back(67);
340  fDecayTo.at(8).push_back(70);
341  fDecayTo.at(8).push_back(71);
342 
343  fBranchingRatios.at(9).push_back(0.16129);
344  fBranchingRatios.at(9).push_back(0.193548);
345  fBranchingRatios.at(9).push_back(0.645161);
346  fDecayTo.at(9).push_back(37);
347  fDecayTo.at(9).push_back(65);
348  fDecayTo.at(9).push_back(72);
349 
350  fBranchingRatios.at(10).push_back(0.245826);
351  fBranchingRatios.at(10).push_back(0.0816327);
352  fBranchingRatios.at(10).push_back(0.20872);
353  fBranchingRatios.at(10).push_back(0.463822);
354  fDecayTo.at(10).push_back(40);
355  fDecayTo.at(10).push_back(56);
356  fDecayTo.at(10).push_back(60);
357  fDecayTo.at(10).push_back(71);
358 
359  fBranchingRatios.at(11).push_back(0.170984);
360  fBranchingRatios.at(11).push_back(0.310881);
361  fBranchingRatios.at(11).push_back(0.518135);
362  fDecayTo.at(11).push_back(29);
363  fDecayTo.at(11).push_back(54);
364  fDecayTo.at(11).push_back(66);
365 
366  fBranchingRatios.at(12).push_back(0.242424);
367  fBranchingRatios.at(12).push_back(0.757576);
368  fDecayTo.at(12).push_back(54);
369  fDecayTo.at(12).push_back(62);
370 
371  fBranchingRatios.at(13).push_back(0.159664);
372  fBranchingRatios.at(13).push_back(0.840336);
373  fDecayTo.at(13).push_back(48);
374  fDecayTo.at(13).push_back(58);
375 
376  fBranchingRatios.at(14).push_back(0.288136);
377  fBranchingRatios.at(14).push_back(0.145763);
378  fBranchingRatios.at(14).push_back(0.118644);
379  fBranchingRatios.at(14).push_back(0.108475);
380  fBranchingRatios.at(14).push_back(0.338983);
381  fDecayTo.at(14).push_back(56);
382  fDecayTo.at(14).push_back(66);
383  fDecayTo.at(14).push_back(70);
384  fDecayTo.at(14).push_back(71);
385  fDecayTo.at(14).push_back(72);
386 
387  fBranchingRatios.at(15).push_back(0.00869263);
388  fBranchingRatios.at(15).push_back(0.087274);
389  fBranchingRatios.at(15).push_back(0.0973574);
390  fBranchingRatios.at(15).push_back(0.288595);
391  fBranchingRatios.at(15).push_back(0.347705);
392  fBranchingRatios.at(15).push_back(0.170376);
393  fDecayTo.at(15).push_back(40);
394  fDecayTo.at(15).push_back(64);
395  fDecayTo.at(15).push_back(65);
396  fDecayTo.at(15).push_back(68);
397  fDecayTo.at(15).push_back(70);
398  fDecayTo.at(15).push_back(71);
399 
400  fBranchingRatios.at(16).push_back(1);
401  fDecayTo.at(16).push_back(65);
402 
403  fBranchingRatios.at(17).push_back(0.341763);
404  fBranchingRatios.at(17).push_back(0.0567327);
405  fBranchingRatios.at(17).push_back(0.164046);
406  fBranchingRatios.at(17).push_back(0.239234);
407  fBranchingRatios.at(17).push_back(0.198223);
408  fDecayTo.at(17).push_back(35);
409  fDecayTo.at(17).push_back(39);
410  fDecayTo.at(17).push_back(51);
411  fDecayTo.at(17).push_back(67);
412  fDecayTo.at(17).push_back(69);
413 
414  fBranchingRatios.at(18).push_back(0.106406);
415  fBranchingRatios.at(18).push_back(0.254103);
416  fBranchingRatios.at(18).push_back(0.0465855);
417  fBranchingRatios.at(18).push_back(0.529381);
418  fBranchingRatios.at(18).push_back(0.0635257);
419  fDecayTo.at(18).push_back(60);
420  fDecayTo.at(18).push_back(61);
421  fDecayTo.at(18).push_back(63);
422  fDecayTo.at(18).push_back(70);
423  fDecayTo.at(18).push_back(72);
424 
425  fBranchingRatios.at(19).push_back(0.0905797);
426  fBranchingRatios.at(19).push_back(0.228261);
427  fBranchingRatios.at(19).push_back(0.181159);
428  fBranchingRatios.at(19).push_back(0.137681);
429  fBranchingRatios.at(19).push_back(0.362319);
430  fDecayTo.at(19).push_back(43);
431  fDecayTo.at(19).push_back(45);
432  fDecayTo.at(19).push_back(52);
433  fDecayTo.at(19).push_back(70);
434  fDecayTo.at(19).push_back(71);
435 
436  fBranchingRatios.at(20).push_back(0.0316414);
437  fBranchingRatios.at(20).push_back(0.0415293);
438  fBranchingRatios.at(20).push_back(0.0362558);
439  fBranchingRatios.at(20).push_back(0.0481213);
440  fBranchingRatios.at(20).push_back(0.0903098);
441  fBranchingRatios.at(20).push_back(0.0929466);
442  fBranchingRatios.at(20).push_back(0.659196);
443  fDecayTo.at(20).push_back(30);
444  fDecayTo.at(20).push_back(32);
445  fDecayTo.at(20).push_back(46);
446  fDecayTo.at(20).push_back(61);
447  fDecayTo.at(20).push_back(64);
448  fDecayTo.at(20).push_back(66);
449  fDecayTo.at(20).push_back(70);
450 
451  fBranchingRatios.at(21).push_back(0.310757);
452  fBranchingRatios.at(21).push_back(0.398406);
453  fBranchingRatios.at(21).push_back(0.290837);
454  fDecayTo.at(21).push_back(64);
455  fDecayTo.at(21).push_back(66);
456  fDecayTo.at(21).push_back(70);
457 
458  fBranchingRatios.at(22).push_back(0.152542);
459  fBranchingRatios.at(22).push_back(0.847458);
460  fDecayTo.at(22).push_back(67);
461  fDecayTo.at(22).push_back(71);
462 
463  fBranchingRatios.at(23).push_back(0.168367);
464  fBranchingRatios.at(23).push_back(0.321429);
465  fBranchingRatios.at(23).push_back(0.510204);
466  fDecayTo.at(23).push_back(52);
467  fDecayTo.at(23).push_back(70);
468  fDecayTo.at(23).push_back(71);
469 
470  fBranchingRatios.at(24).push_back(0.0276437);
471  fBranchingRatios.at(24).push_back(0.078982);
472  fBranchingRatios.at(24).push_back(0.0245722);
473  fBranchingRatios.at(24).push_back(0.162352);
474  fBranchingRatios.at(24).push_back(0.184291);
475  fBranchingRatios.at(24).push_back(0.438789);
476  fBranchingRatios.at(24).push_back(0.0833699);
477  fDecayTo.at(24).push_back(36);
478  fDecayTo.at(24).push_back(53);
479  fDecayTo.at(24).push_back(62);
480  fDecayTo.at(24).push_back(64);
481  fDecayTo.at(24).push_back(70);
482  fDecayTo.at(24).push_back(71);
483  fDecayTo.at(24).push_back(72);
484 
485  fBranchingRatios.at(25).push_back(0.0163934);
486  fBranchingRatios.at(25).push_back(0.336276);
487  fBranchingRatios.at(25).push_back(0.226986);
488  fBranchingRatios.at(25).push_back(0.420345);
489  fDecayTo.at(25).push_back(43);
490  fDecayTo.at(25).push_back(67);
491  fDecayTo.at(25).push_back(68);
492  fDecayTo.at(25).push_back(70);
493 
494  fBranchingRatios.at(26).push_back(0.184332);
495  fBranchingRatios.at(26).push_back(0.0460829);
496  fBranchingRatios.at(26).push_back(0.460829);
497  fBranchingRatios.at(26).push_back(0.078341);
498  fBranchingRatios.at(26).push_back(0.230415);
499  fDecayTo.at(26).push_back(53);
500  fDecayTo.at(26).push_back(54);
501  fDecayTo.at(26).push_back(60);
502  fDecayTo.at(26).push_back(61);
503  fDecayTo.at(26).push_back(71);
504 
505  fBranchingRatios.at(27).push_back(0.0147145);
506  fBranchingRatios.at(27).push_back(0.0170689);
507  fBranchingRatios.at(27).push_back(0.0500294);
508  fBranchingRatios.at(27).push_back(0.329606);
509  fBranchingRatios.at(27).push_back(0.588582);
510  fDecayTo.at(27).push_back(36);
511  fDecayTo.at(27).push_back(46);
512  fDecayTo.at(27).push_back(56);
513  fDecayTo.at(27).push_back(67);
514  fDecayTo.at(27).push_back(68);
515 
516  fBranchingRatios.at(28).push_back(1);
517  fDecayTo.at(28).push_back(70);
518 
519  fBranchingRatios.at(29).push_back(0.280702);
520  fBranchingRatios.at(29).push_back(0.0935673);
521  fBranchingRatios.at(29).push_back(0.0409357);
522  fBranchingRatios.at(29).push_back(0.584795);
523  fDecayTo.at(29).push_back(63);
524  fDecayTo.at(29).push_back(66);
525  fDecayTo.at(29).push_back(68);
526  fDecayTo.at(29).push_back(70);
527 
528  fBranchingRatios.at(30).push_back(0.393701);
529  fBranchingRatios.at(30).push_back(0.283465);
530  fBranchingRatios.at(30).push_back(0.188976);
531  fBranchingRatios.at(30).push_back(0.133858);
532  fDecayTo.at(30).push_back(61);
533  fDecayTo.at(30).push_back(67);
534  fDecayTo.at(30).push_back(71);
535  fDecayTo.at(30).push_back(72);
536 
537  fBranchingRatios.at(31).push_back(0.0477737);
538  fBranchingRatios.at(31).push_back(0.194805);
539  fBranchingRatios.at(31).push_back(0.0245826);
540  fBranchingRatios.at(31).push_back(0.269017);
541  fBranchingRatios.at(31).push_back(0.463822);
542  fDecayTo.at(31).push_back(45);
543  fDecayTo.at(31).push_back(60);
544  fDecayTo.at(31).push_back(65);
545  fDecayTo.at(31).push_back(71);
546  fDecayTo.at(31).push_back(72);
547 
548  fBranchingRatios.at(32).push_back(0.105769);
549  fBranchingRatios.at(32).push_back(0.129808);
550  fBranchingRatios.at(32).push_back(0.0528846);
551  fBranchingRatios.at(32).push_back(0.480769);
552  fBranchingRatios.at(32).push_back(0.230769);
553  fDecayTo.at(32).push_back(46);
554  fDecayTo.at(32).push_back(56);
555  fDecayTo.at(32).push_back(66);
556  fDecayTo.at(32).push_back(70);
557  fDecayTo.at(32).push_back(71);
558 
559  fBranchingRatios.at(33).push_back(0.21875);
560  fBranchingRatios.at(33).push_back(0.78125);
561  fDecayTo.at(33).push_back(67);
562  fDecayTo.at(33).push_back(71);
563 
564  fBranchingRatios.at(34).push_back(0.102011);
565  fBranchingRatios.at(34).push_back(0.179598);
566  fBranchingRatios.at(34).push_back(0.718391);
567  fDecayTo.at(34).push_back(43);
568  fDecayTo.at(34).push_back(61);
569  fDecayTo.at(34).push_back(66);
570 
571  fBranchingRatios.at(35).push_back(0.00826446);
572  fBranchingRatios.at(35).push_back(0.393546);
573  fBranchingRatios.at(35).push_back(0.334514);
574  fBranchingRatios.at(35).push_back(0.263676);
575  fDecayTo.at(35).push_back(64);
576  fDecayTo.at(35).push_back(67);
577  fDecayTo.at(35).push_back(68);
578  fDecayTo.at(35).push_back(70);
579 
580  fBranchingRatios.at(36).push_back(0.056338);
581  fBranchingRatios.at(36).push_back(0.704225);
582  fBranchingRatios.at(36).push_back(0.239437);
583  fDecayTo.at(36).push_back(51);
584  fDecayTo.at(36).push_back(70);
585  fDecayTo.at(36).push_back(71);
586 
587  fBranchingRatios.at(37).push_back(0.21875);
588  fBranchingRatios.at(37).push_back(0.78125);
589  fDecayTo.at(37).push_back(67);
590  fDecayTo.at(37).push_back(70);
591 
592  fBranchingRatios.at(38).push_back(0.181818);
593  fBranchingRatios.at(38).push_back(0.757576);
594  fBranchingRatios.at(38).push_back(0.0606061);
595  fDecayTo.at(38).push_back(66);
596  fDecayTo.at(38).push_back(71);
597  fDecayTo.at(38).push_back(72);
598 
599  fBranchingRatios.at(39).push_back(0.157258);
600  fBranchingRatios.at(39).push_back(0.403226);
601  fBranchingRatios.at(39).push_back(0.237903);
602  fBranchingRatios.at(39).push_back(0.201613);
603  fDecayTo.at(39).push_back(62);
604  fDecayTo.at(39).push_back(70);
605  fDecayTo.at(39).push_back(71);
606  fDecayTo.at(39).push_back(72);
607 
608  fBranchingRatios.at(40).push_back(0.0740741);
609  fBranchingRatios.at(40).push_back(0.925926);
610  fDecayTo.at(40).push_back(52);
611  fDecayTo.at(40).push_back(72);
612 
613  fBranchingRatios.at(41).push_back(0.0535714);
614  fBranchingRatios.at(41).push_back(0.35119);
615  fBranchingRatios.at(41).push_back(0.595238);
616  fDecayTo.at(41).push_back(67);
617  fDecayTo.at(41).push_back(68);
618  fDecayTo.at(41).push_back(70);
619 
620  fBranchingRatios.at(42).push_back(0.00816803);
621  fBranchingRatios.at(42).push_back(0.0583431);
622  fBranchingRatios.at(42).push_back(0.350058);
623  fBranchingRatios.at(42).push_back(0.583431);
624  fDecayTo.at(42).push_back(49);
625  fDecayTo.at(42).push_back(62);
626  fDecayTo.at(42).push_back(71);
627  fDecayTo.at(42).push_back(72);
628 
629  fBranchingRatios.at(43).push_back(0.0961538);
630  fBranchingRatios.at(43).push_back(0.423077);
631  fBranchingRatios.at(43).push_back(0.480769);
632  fDecayTo.at(43).push_back(66);
633  fDecayTo.at(43).push_back(67);
634  fDecayTo.at(43).push_back(68);
635 
636  fBranchingRatios.at(44).push_back(0.450549);
637  fBranchingRatios.at(44).push_back(0.549451);
638  fDecayTo.at(44).push_back(69);
639  fDecayTo.at(44).push_back(72);
640 
641  fBranchingRatios.at(45).push_back(0.469925);
642  fBranchingRatios.at(45).push_back(0.0836466);
643  fBranchingRatios.at(45).push_back(0.446429);
644  fDecayTo.at(45).push_back(61);
645  fDecayTo.at(45).push_back(65);
646  fDecayTo.at(45).push_back(72);
647 
648  fBranchingRatios.at(46).push_back(0.0408163);
649  fBranchingRatios.at(46).push_back(0.510204);
650  fBranchingRatios.at(46).push_back(0.44898);
651  fDecayTo.at(46).push_back(67);
652  fDecayTo.at(46).push_back(70);
653  fDecayTo.at(46).push_back(71);
654 
655  fBranchingRatios.at(47).push_back(1);
656  fDecayTo.at(47).push_back(72);
657 
658  fBranchingRatios.at(48).push_back(0.641026);
659  fBranchingRatios.at(48).push_back(0.358974);
660  fDecayTo.at(48).push_back(58);
661  fDecayTo.at(48).push_back(69);
662 
663  fBranchingRatios.at(49).push_back(0.0458015);
664  fBranchingRatios.at(49).push_back(0.954198);
665  fDecayTo.at(49).push_back(66);
666  fDecayTo.at(49).push_back(70);
667 
668  fBranchingRatios.at(50).push_back(0.401639);
669  fBranchingRatios.at(50).push_back(0.188525);
670  fBranchingRatios.at(50).push_back(0.409836);
671  fDecayTo.at(50).push_back(61);
672  fDecayTo.at(50).push_back(69);
673  fDecayTo.at(50).push_back(72);
674 
675  fBranchingRatios.at(51).push_back(0.0188679);
676  fBranchingRatios.at(51).push_back(0.173585);
677  fBranchingRatios.at(51).push_back(0.754717);
678  fBranchingRatios.at(51).push_back(0.0528302);
679  fDecayTo.at(51).push_back(61);
680  fDecayTo.at(51).push_back(67);
681  fDecayTo.at(51).push_back(71);
682  fDecayTo.at(51).push_back(72);
683 
684  fBranchingRatios.at(52).push_back(0.0100849);
685  fBranchingRatios.at(52).push_back(0.00796178);
686  fBranchingRatios.at(52).push_back(0.530786);
687  fBranchingRatios.at(52).push_back(0.451168);
688  fDecayTo.at(52).push_back(59);
689  fDecayTo.at(52).push_back(68);
690  fDecayTo.at(52).push_back(70);
691  fDecayTo.at(52).push_back(71);
692 
693  fBranchingRatios.at(53).push_back(0.0379902);
694  fBranchingRatios.at(53).push_back(0.0490196);
695  fBranchingRatios.at(53).push_back(0.612745);
696  fBranchingRatios.at(53).push_back(0.300245);
697  fDecayTo.at(53).push_back(67);
698  fDecayTo.at(53).push_back(70);
699  fDecayTo.at(53).push_back(71);
700  fDecayTo.at(53).push_back(72);
701 
702  fBranchingRatios.at(54).push_back(0.106557);
703  fBranchingRatios.at(54).push_back(0.819672);
704  fBranchingRatios.at(54).push_back(0.0737705);
705  fDecayTo.at(54).push_back(59);
706  fDecayTo.at(54).push_back(68);
707  fDecayTo.at(54).push_back(70);
708 
709  fBranchingRatios.at(55).push_back(0.699301);
710  fBranchingRatios.at(55).push_back(0.300699);
711  fDecayTo.at(55).push_back(64);
712  fDecayTo.at(55).push_back(70);
713 
714  fBranchingRatios.at(56).push_back(1);
715  fDecayTo.at(56).push_back(71);
716 
717  fBranchingRatios.at(57).push_back(1);
718  fDecayTo.at(57).push_back(72);
719 
720  fBranchingRatios.at(58).push_back(0.888099);
721  fBranchingRatios.at(58).push_back(0.111901);
722  fDecayTo.at(58).push_back(69);
723  fDecayTo.at(58).push_back(72);
724 
725  fBranchingRatios.at(59).push_back(0.00647298);
726  fBranchingRatios.at(59).push_back(0.752672);
727  fBranchingRatios.at(59).push_back(0.165588);
728  fBranchingRatios.at(59).push_back(0.0752672);
729  fDecayTo.at(59).push_back(65);
730  fDecayTo.at(59).push_back(70);
731  fDecayTo.at(59).push_back(71);
732  fDecayTo.at(59).push_back(72);
733 
734  fBranchingRatios.at(60).push_back(0.0708556);
735  fBranchingRatios.at(60).push_back(0.668449);
736  fBranchingRatios.at(60).push_back(0.260695);
737  fDecayTo.at(60).push_back(65);
738  fDecayTo.at(60).push_back(71);
739  fDecayTo.at(60).push_back(72);
740 
741  fBranchingRatios.at(61).push_back(0.166667);
742  fBranchingRatios.at(61).push_back(0.833333);
743  fDecayTo.at(61).push_back(69);
744  fDecayTo.at(61).push_back(72);
745 
746  fBranchingRatios.at(62).push_back(0.0898551);
747  fBranchingRatios.at(62).push_back(0.57971);
748  fBranchingRatios.at(62).push_back(0.330435);
749  fDecayTo.at(62).push_back(67);
750  fDecayTo.at(62).push_back(68);
751  fDecayTo.at(62).push_back(70);
752 
753  fBranchingRatios.at(63).push_back(0.813008);
754  fBranchingRatios.at(63).push_back(0.186992);
755  fDecayTo.at(63).push_back(71);
756  fDecayTo.at(63).push_back(72);
757 
758  fBranchingRatios.at(64).push_back(0.29078);
759  fBranchingRatios.at(64).push_back(0.70922);
760  fDecayTo.at(64).push_back(70);
761  fDecayTo.at(64).push_back(71);
762 
763  fBranchingRatios.at(65).push_back(0.05);
764  fBranchingRatios.at(65).push_back(0.08);
765  fBranchingRatios.at(65).push_back(0.5);
766  fBranchingRatios.at(65).push_back(0.37);
767  fDecayTo.at(65).push_back(69);
768  fDecayTo.at(65).push_back(70);
769  fDecayTo.at(65).push_back(71);
770  fDecayTo.at(65).push_back(72);
771 
772  fBranchingRatios.at(66).push_back(0.398406);
773  fBranchingRatios.at(66).push_back(0.310757);
774  fBranchingRatios.at(66).push_back(0.290837);
775  fDecayTo.at(66).push_back(70);
776  fDecayTo.at(66).push_back(71);
777  fDecayTo.at(66).push_back(72);
778 
779  fBranchingRatios.at(67).push_back(0.819672);
780  fBranchingRatios.at(67).push_back(0.180328);
781  fDecayTo.at(67).push_back(70);
782  fDecayTo.at(67).push_back(71);
783 
784  fBranchingRatios.at(68).push_back(0.186992);
785  fBranchingRatios.at(68).push_back(0.813008);
786  fDecayTo.at(68).push_back(70);
787  fDecayTo.at(68).push_back(71);
788 
789  fBranchingRatios.at(69).push_back(1);
790  fDecayTo.at(69).push_back(72);
791 
792  fBranchingRatios.at(70).push_back(1);
793  fDecayTo.at(70).push_back(71);
794 
795  fBranchingRatios.at(71).push_back(1);
796  fDecayTo.at(71).push_back(72);
797 
798  // Info from table VII (http://prc.aps.org/pdf/PRC/v58/i6/p3677_1)
799  // in MeV
800  double startEnergyLevels[] = { 2.281, 2.752, 2.937,
801  3.143, 3.334, 3.569, 3.652, 3.786, 3.861, 4.067, 4.111, 4.267,
802  4.364, 4.522, 4.655, 4.825, 5.017, 5.080, 5.223, 5.696, 6.006 };
803 
804  fStartEnergyLevels.insert(fStartEnergyLevels.end(), &startEnergyLevels[0],
805  &startEnergyLevels[fNumberOfStartLevels]);
806 
807  double b[] = { 0.9, 1.5, 0.11, 0.06, 0.04, 0.01,
808  0.16, 0.26, 0.01, 0.05, 0.11, 0.29,
809  3.84, 0.31, 0.38, 0.47, 0.36, 0.23,
810  0.03, 0.11, 0.13 };
811 
812  fB.insert(fB.end(), &b[0], &b[fNumberOfStartLevels]);
813 
814  double energyLevels[] = { 7.4724, 6.2270, 5.06347,
815  4.99294, 4.8756, 4.87255, 4.78865, 4.744093, 4.53706,
816  4.47299, 4.41936, 4.39588, 4.3840, 4.3656, 4.28052,
817  4.25362, 4.21307, 4.18003, 4.14901, 4.11084, 4.10446,
818  4.02035, 3.92390, 3.88792, 3.86866, 3.840228, 3.82143,
819  3.79757, 3.76779, 3.73848, 3.663739, 3.62995, 3.59924,
820  3.55697, 3.48621, 3.439144, 3.41434, 3.39363, 3.36803,
821  3.22867, 3.15381, 3.14644, 3.12836, 3.109721, 3.1002,
822  3.02795, 2.98587, 2.9508, 2.87901, 2.80788, 2.7874,
823  2.786644, 2.75672, 2.74691, 2.730372, 2.625990, 2.57593,
824  2.558, 2.54277, 2.419171, 2.397165, 2.290493, 2.289871,
825  2.26040, 2.103668, 2.069809, 2.047354, 1.959068, 1.643639,
826  0.891398, 0.8001427, 0.0298299, 0.0 };
827 
828  fEnergyLevels.insert(fEnergyLevels.end(), &energyLevels[0],
829  &energyLevels[fNumberOfLevels]);
830 
831  // I have to put checks that the arrays really
832  // have as many elements as they should
833 
834  }
evgen::NueAr40CCGenerator::fStartEnergyLevels
std::vector< double > fStartEnergyLevels
Definition: NueAr40CCGenerator.h:86
evgen::NueAr40CCGenerator::fEnergyLevels
std::vector< double > fEnergyLevels
Definition: NueAr40CCGenerator.h:88
evgen::NueAr40CCGenerator::fBranchingRatios
std::vector< std::vector< double > > fBranchingRatios
Definition: NueAr40CCGenerator.h:84
evgen::NueAr40CCGenerator::fB
std::vector< double > fB
Definition: NueAr40CCGenerator.h:87
evgen::NueAr40CCGenerator::fDecayTo
std::vector< std::vector< int > > fDecayTo
Definition: NueAr40CCGenerator.h:85
bool evgen::NueAr40CCGenerator::ProcessOneNeutrino ( simb::MCTruth &  truth,
double  neutrinoEnergy,
double  neutrinoTime,
CLHEP::HepRandomEngine &  engine 
) const
private

Definition at line 854 of file NueAr40CCGenerator.cxx.

References simb::MCTruth::Add(), simb::MCParticle::AddTrajectoryPoint(), CalculateCrossSections(), fBranchingRatios, fDecayTo, fEnergyLevels, fNumberOfLevels, fStartEnergyLevels, GetIsotropicDirection(), GetUniformPosition(), simb::kCC, simb::kCCQE, simb::kQE, n, simb::MCTruth::NParticles(), and simb::MCTruth::SetNeutrino().

Referenced by CreateKinematicsVector().

857  {
858 
859  CLHEP::RandFlat randFlat(engine);
860 
861  int highestLevel = 0;
862  std::vector< double > levelCrossSections =
863  CalculateCrossSections(neutrinoEnergy, highestLevel);
864 
865  double totalCrossSection = 0;
866  // Calculating total cross section
867  for (std::vector< double >::iterator crossSection =
868  levelCrossSections.begin();
869  crossSection != levelCrossSections.end(); ++crossSection)
870  totalCrossSection += *crossSection;
871 
872  if (totalCrossSection == 0)
873  return false;
874 
875  std::vector< double > startLevelProbabilities;
876  // Calculating each level's probability
877  for (std::vector< double >::iterator crossSection =
878  levelCrossSections.begin();
879  crossSection != levelCrossSections.end(); ++crossSection)
880  startLevelProbabilities.push_back((*crossSection)/totalCrossSection);
881 
882  double randomNumber = randFlat.fire();
883  double tprob = 0;
884  int chosenStartLevel = -1;
885  // Picking a starting level
886  for (int level = 0; level < highestLevel; ++level)
887  {
888  if (randomNumber < (startLevelProbabilities.at(level) + tprob))
889  {
890  chosenStartLevel = level;
891  break;
892  }
893  tprob += startLevelProbabilities.at(level);
894  }
895 
896  int lastLevel = -1;
897  int level = -1;
898 
899  // Time delays
900  // Seems like this is not implemented yet
901  std::vector< double > levelDelay;
902  for (int n = 0; n < fNumberOfLevels; ++n) levelDelay.push_back(0.0);
903 
904  // The highest n for which fStartEnergyLevels.at(chosenStartLevel)
905  // is higher than fEnergyLevels.at(n)
906  int highestHigher = 0;
907  // The lowest n for which fStartEnergyLevels.at(chosenStartLevel)
908  // is lower than fEnergyLevels.at(n)
909  int lowestLower = 0;
910 
911  // Finding lowestLower and highestHigher
912  for (int n = 0; n < fNumberOfLevels; ++n)
913  {
914  if (fStartEnergyLevels.at(chosenStartLevel) < fEnergyLevels.at(n))
915  lowestLower = n;
916  if (fStartEnergyLevels.at(chosenStartLevel) > fEnergyLevels.at(n))
917  {
918  highestHigher = n;
919  break;
920  }
921  }
922 
923  if (std::abs(fStartEnergyLevels.at(chosenStartLevel) -
924  fEnergyLevels.at(lowestLower))
925  < std::abs(fStartEnergyLevels.at(chosenStartLevel)
926  - fEnergyLevels.at(highestHigher)))
927  {
928  lastLevel = lowestLower;
929  level = lowestLower;
930  }
931  // If the chosen start level energy is closest to the lowest level energy
932  // that it's lower than than the highest level energy
933  // that it's higher than, it starts at the level of
934  // the lowest level energy that it's lower than
935 
936  if (std::abs(fStartEnergyLevels.at(chosenStartLevel)
937  - fEnergyLevels.at(highestHigher))
938  < std::abs(fStartEnergyLevels.at(chosenStartLevel)
939  - fEnergyLevels.at(lowestLower)))
940  {
941  lastLevel = highestHigher;
942  level = highestHigher;
943  }
944  // If the chosen start level energy is closest to the highest level energy
945  // that it's higher than than the lowest level energy that it's lower than,
946  // it starts at the level of the highest level energy that it's higher than
947 
948  std::vector< double > vertex = GetUniformPosition(engine);
949 
950  std::vector< double > neutrinoDirection = GetIsotropicDirection(engine);
951 
952 
953 
954  //
955  // Add the first particle (the neutrino) to the truth list...
956  //
957 
958  // NOTE: all of the values below calculated for the neutrino should be checked!
959 
960  // Primary particles have negative IDs
961  int neutrinoTrackID = -1*(truth.NParticles() + 1);
962  std::string primary("primary");
963 
964  int nuePDG = 12;
965  double neutrinoP = neutrinoEnergy/1000.0; // use GeV...
966  double neutrinoPx = neutrinoDirection.at(0)*neutrinoP;
967  double neutrinoPy = neutrinoDirection.at(1)*neutrinoP;
968  double neutrinoPz = neutrinoDirection.at(2)*neutrinoP;
969 
970 
971 
972  // Create the neutrino:
973  // * set the mother to -1 since it is primary
974  // * set the mass to something < 0 so that the constructor looks up the mass from the pdg tables
975  // * set the status bit to 0 so that geant doesn't waste any CPU tracking the neutrino
976  simb::MCParticle neutrino(neutrinoTrackID, nuePDG, primary, -1, -1, 0);
977 
978  TLorentzVector neutrinoPosition(vertex.at(0), vertex.at(1),
979  vertex.at(2), neutrinoTime);
980  TLorentzVector neutrinoMomentum(neutrinoPx, neutrinoPy,
981  neutrinoPz, neutrinoEnergy/1000.0);
982  neutrino.AddTrajectoryPoint(neutrinoPosition, neutrinoMomentum);
983 
984  truth.Add(neutrino);
985 
986 
987 
988  //
989  // Adding the electron to truth
990  //
991 
992  // In MeV
993  double electronEnergy = neutrinoEnergy -
994  (fEnergyLevels.at(lastLevel) + 1.5);
995  double electronEnergyGeV = electronEnergy/1000.0;
996 
997  double electronM = 0.000511;
998  double electronP = std::sqrt(std::pow(electronEnergyGeV,2)
999  - std::pow(electronM,2));
1000 
1001  // For the moment, choose an isotropic direction for the electron.
1002  std::vector< double > electronDirection = GetIsotropicDirection(engine);
1003  double electronPx = electronDirection.at(0)*electronP;
1004  double electronPy = electronDirection.at(1)*electronP;
1005  double electronPz = electronDirection.at(2)*electronP;
1006 
1007  // Primary particles have negative IDs
1008  int trackID = -1*(truth.NParticles() + 1);
1009  int electronPDG = 11;
1010  simb::MCParticle electron(trackID, electronPDG, primary);
1011 
1012  TLorentzVector electronPosition(vertex.at(0), vertex.at(1),
1013  vertex.at(2), neutrinoTime);
1014  TLorentzVector electronMomentum(electronPx, electronPy,
1015  electronPz, electronEnergyGeV);
1016  electron.AddTrajectoryPoint(electronPosition, electronMomentum);
1017 
1018  truth.Add(electron);
1019 
1020  double ttime = neutrinoTime;
1021  int noMoreDecay = 0;
1022 
1023  // Level loop
1024  int groundLevel = fNumberOfLevels - 1;
1025  while (level != groundLevel)
1026  {
1027 
1028  double rl = randFlat.fire();
1029 
1030  int decayNum = 0;
1031 
1032  tprob = 0; // Used this variable above for cross section
1033 
1034  // Decay loop
1035  for (unsigned int iLevel = 0;
1036  iLevel < fBranchingRatios.at(level).size(); ++iLevel)
1037  {
1038 
1039  if (rl < (fBranchingRatios.at(level).at(iLevel) + tprob))
1040  {
1041 
1042  // We have a decay
1043 
1044  level = fDecayTo.at(level).at(decayNum);
1045 
1046  // Really should be using a map
1047  double gammaEnergy = fEnergyLevels.at(lastLevel) -
1048  fEnergyLevels.at(level);
1049  double gammaEnergyGeV = gammaEnergy/1000;
1050 
1051  std::vector< double > gammaDirection = GetIsotropicDirection(engine);
1052  double gammaPx = gammaDirection.at(0)*gammaEnergyGeV;
1053  double gammaPy = gammaDirection.at(1)*gammaEnergyGeV;
1054  double gammaPz = gammaDirection.at(2)*gammaEnergyGeV;
1055 
1056  //double gammaM = 0.0; // unused
1057 
1058  double gammaTime = (-TMath::Log(randFlat.fire())/
1059  (1/(levelDelay.at(lastLevel)))) + ttime;
1060 
1061  // Adding the gamma to truth
1062 
1063  // Primary particles have negative IDs
1064  trackID = -1*(truth.NParticles() + 1);
1065  int gammaPDG = 22;
1066  simb::MCParticle gamma(trackID, gammaPDG, primary); // NOTE: should the gammas be primary?
1067 
1068  TLorentzVector gammaPosition(vertex.at(0), vertex.at(1),
1069  vertex.at(2), neutrinoTime);
1070  TLorentzVector gammaMomentum(gammaPx, gammaPy,
1071  gammaPz, gammaEnergyGeV);
1072  gamma.AddTrajectoryPoint(gammaPosition, gammaMomentum);
1073 
1074  truth.Add(gamma);
1075 
1076  lastLevel = level;
1077  ttime = gammaTime;
1078 
1079  break;
1080 
1081  }
1082 
1083  if ((tprob + fBranchingRatios.at(level).at(iLevel)) > 1) {
1084  std::cout << "(tprob + *iLevel) > 1" << std::endl;
1085  noMoreDecay = 1; // If it doesn't do any more gamma decay
1086  break;
1087  }
1088 
1089  ++decayNum;
1090  tprob += fBranchingRatios.at(level).at(iLevel);
1091 
1092  } // End of decay loop
1093 
1094  if (noMoreDecay == 1) break;
1095 
1096  } // End of level loop
1097 
1098  if (level != 72)
1099  {
1100  std::cout << "level != 72" << std::endl;
1101  std::cout << "level = " << level << std::endl;
1102  }
1103 
1104  // Set the neutrino for the MCTruth object:
1105  // NOTE: currently these parameters are all pretty much a guess...
1106  truth.SetNeutrino(simb::kCC,
1107  simb::kQE, // not sure what the mode should be here, assumed that these are all QE???
1108  simb::kCCQE, // not sure what the int_type should be either...
1109  0, // target is AR40? not sure how to specify that...
1110  0, // nucleon pdg ???
1111  0, // quark pdg ???
1112  -1.0, // W ??? - not sure how to calculate this from the above
1113  -1.0, // X ??? - not sure how to calculate this from the above
1114  -1.0, // Y ??? - not sure how to calculate this from the above
1115  -1.0); // Qsqr ??? - not sure how to calculate this from the above
1116 
1117  return true;
1118 
1119  }
evgen::NueAr40CCGenerator::fStartEnergyLevels
std::vector< double > fStartEnergyLevels
Definition: NueAr40CCGenerator.h:86
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:24
evgen::NueAr40CCGenerator::GetUniformPosition
std::vector< double > GetUniformPosition(CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:120
evgen::NueAr40CCGenerator::fEnergyLevels
std::vector< double > fEnergyLevels
Definition: NueAr40CCGenerator.h:88
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:72
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:100
evgen::NueAr40CCGenerator::fBranchingRatios
std::vector< std::vector< double > > fBranchingRatios
Definition: NueAr40CCGenerator.h:84
simb::MCTruth::Add
void Add(simb::MCParticle &part)
Definition: MCTruth.h:77
evgen::NueAr40CCGenerator::GetIsotropicDirection
std::vector< double > GetIsotropicDirection(CLHEP::HepRandomEngine &engine) const
Definition: NueAr40CCGenerator.cxx:97
evgen::NueAr40CCGenerator::CalculateCrossSections
std::vector< double > CalculateCrossSections(double neutrinoEnergy, int &highestLevel) const
Definition: NueAr40CCGenerator.cxx:1127
n
Char_t n[5]
Definition: comparison_ascii.C:55
evgen::NueAr40CCGenerator::fDecayTo
std::vector< std::vector< int > > fDecayTo
Definition: NueAr40CCGenerator.h:85
vertex
vertex reconstruction
Definition: AggregateVertex_module.cc:43
void evgen::NueAr40CCGenerator::ReadNeutrinoSpectrum ( )
private

Definition at line 215 of file NueAr40CCGenerator.cxx.

References fEnergyProbabilityMap, and fEnergySpectrumFileName.

Referenced by NueAr40CCGenerator().

216  {
217 
218  cet::search_path searchPath("FW_SEARCH_PATH");
219  std::string directoryName = "SupernovaNeutrinos/" +
220  fEnergySpectrumFileName;
221 
222  std::string fullName;
223  searchPath.find_file(directoryName, fullName);
224 
225  if (fullName.empty())
226  throw cet::exception("NueAr40CCGenerator")
227  << "Cannot find file with neutrino energy spectrum "
228  << fullName << '\n';
229 
230  TFile energySpectrumFile(fullName.c_str(), "READ");
231 
232  std::string energySpectrumGraphName = "NueSpectrum";
233  TGraph *energySpectrumGraph =
234  dynamic_cast< TGraph* >(energySpectrumFile.Get
235  (energySpectrumGraphName.c_str()));
236 
237  double integral = 0.0;
238  int numberOfPoints = energySpectrumGraph->GetN();
239  double *energyValues = energySpectrumGraph->GetX();
240  double *fluxValues = energySpectrumGraph->GetY();
241  for (int point = 0; point < numberOfPoints; ++point)
242  integral += fluxValues[point];
243 
244  double probability = 0.0;
245  for (int point = 0; point < numberOfPoints; ++point)
246  {
247  probability += fluxValues[point]/integral;
248  fEnergyProbabilityMap.insert(std::make_pair(energyValues[point],
249  probability));
250  }
251 
252  }
evgen::NueAr40CCGenerator::fEnergyProbabilityMap
std::map< double, double > fEnergyProbabilityMap
Definition: NueAr40CCGenerator.h:78
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 115 of file NueAr40CCGenerator.h.

Referenced by GetUniformPosition(), and NueAr40CCGenerator().

bool evgen::NueAr40CCGenerator::fAllowZeroNeutrinos
private

Definition at line 104 of file NueAr40CCGenerator.h.

Referenced by GetNumberOfNeutrinos().

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

Definition at line 87 of file NueAr40CCGenerator.h.

Referenced by CalculateCrossSections(), and InitializeVectors().

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

Definition at line 84 of file NueAr40CCGenerator.h.

Referenced by InitializeVectors(), and ProcessOneNeutrino().

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

Definition at line 85 of file NueAr40CCGenerator.h.

Referenced by InitializeVectors(), and ProcessOneNeutrino().

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

Definition at line 88 of file NueAr40CCGenerator.h.

Referenced by InitializeVectors(), and ProcessOneNeutrino().

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

Definition at line 78 of file NueAr40CCGenerator.h.

Referenced by GetNeutrinoEnergy(), and ReadNeutrinoSpectrum().

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

Definition at line 96 of file NueAr40CCGenerator.h.

Referenced by ReadNeutrinoSpectrum().

bool evgen::NueAr40CCGenerator::fMonoenergeticNeutrinos
private

Definition at line 91 of file NueAr40CCGenerator.h.

Referenced by GetNeutrinoEnergy(), and NueAr40CCGenerator().

double evgen::NueAr40CCGenerator::fNeutrinoEnergy
private

Definition at line 93 of file NueAr40CCGenerator.h.

Referenced by GetNeutrinoEnergy().

double evgen::NueAr40CCGenerator::fNeutrinoTimeBegin
private

Definition at line 110 of file NueAr40CCGenerator.h.

Referenced by GetNeutrinoTime().

double evgen::NueAr40CCGenerator::fNeutrinoTimeEnd
private

Definition at line 111 of file NueAr40CCGenerator.h.

Referenced by GetNeutrinoTime().

int evgen::NueAr40CCGenerator::fNumberOfLevels
private

Definition at line 80 of file NueAr40CCGenerator.h.

Referenced by InitializeVectors(), and ProcessOneNeutrino().

int evgen::NueAr40CCGenerator::fNumberOfNeutrinos
private

Definition at line 107 of file NueAr40CCGenerator.h.

Referenced by GetNumberOfNeutrinos().

int evgen::NueAr40CCGenerator::fNumberOfStartLevels
private

Definition at line 81 of file NueAr40CCGenerator.h.

Referenced by CalculateCrossSections(), and InitializeVectors().

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

Definition at line 86 of file NueAr40CCGenerator.h.

Referenced by CalculateCrossSections(), InitializeVectors(), and ProcessOneNeutrino().

bool evgen::NueAr40CCGenerator::fUsePoissonDistribution
private

Definition at line 100 of file NueAr40CCGenerator.h.

Referenced by GetNumberOfNeutrinos().

The documentation for this class was generated from the following files: