#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 32 of file NueAr40CCGenerator.h.

Constructor & Destructor Documentation

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

Definition at line 37 of file NueAr40CCGenerator.cxx.

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

     : fNumberOfLevels(73)
     , fNumberOfStartLevels(21)
     , fBranchingRatios(fNumberOfLevels)
     , fDecayTo(fNumberOfLevels)
     , fMonoenergeticNeutrinos(parameterSet.get<bool>("MonoenergeticNeutrinos"))
     , fNeutrinoEnergy(parameterSet.get<double>("NeutrinoEnergy"))
     , fEnergySpectrumFileName(parameterSet.get<std::string>("EnergySpectrumFileName"))
     , fUsePoissonDistribution(parameterSet.get<bool>("UsePoissonDistribution"))
     , fAllowZeroNeutrinos(parameterSet.get<bool>("AllowZeroNeutrinos"))
     , fNumberOfNeutrinos(parameterSet.get<int>("NumberOfNeutrinos"))
     , fNeutrinoTimeBegin(parameterSet.get<double>("NeutrinoTimeBegin"))
     , fNeutrinoTimeEnd(parameterSet.get<double>("NeutrinoTimeEnd"))
   {
 
     fActiveVolume.push_back(parameterSet.get<std::vector<double>>("ActiveVolume0"));
     fActiveVolume.push_back(parameterSet.get<std::vector<double>>("ActiveVolume1"));
 
     if (!fMonoenergeticNeutrinos) ReadNeutrinoSpectrum();
 
     InitializeVectors();
   }

Member Function Documentation

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

private

Definition at line 1053 of file NueAr40CCGenerator.cxx.

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

Referenced by ProcessOneNeutrino().

   {
 
     highestLevel = 0;
     std::vector<double> levelCrossSections;
 
     // Loop through energy levels, if neutrino has enough energy,
     // calculate cross section
     for (int level = 0; level < fNumberOfStartLevels; ++level) {
       // Electron energy in keV
       double w = (neutrinoEnergy - (fStartEnergyLevels.at(level) + 1.5)) * 1000;
 
       double sigma = 0.0;
       if (neutrinoEnergy > (fStartEnergyLevels.at(level) + 1.5) && w >= 511.) {
         ++highestLevel;
         for (int n = 0; n <= level; ++n) {
           // Electron energy for level n to use in cross section
           w = (neutrinoEnergy - (fStartEnergyLevels.at(n) + 1.5)) * 1000;
           // Electron momentum in keV/c
           double p = std::sqrt(pow(w, 2) - pow(511.0, 2));
           // Fermi function approximation
           double f = std::sqrt(3.0634 + (0.6814 / (w - 1)));
           // In cm^2*10^-42
           sigma += 1.6e-8 * (p * w * f * fB.at(n));
         }
       }
       levelCrossSections.push_back(sigma);
     }
 
     return levelCrossSections;
   }

void evgen::NueAr40CCGenerator::CreateKinematicsVector	(	simb::MCTruth &	truth,
		CLHEP::HepRandomEngine &	engine
	)		const

private

Definition at line 798 of file NueAr40CCGenerator.cxx.

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

Referenced by Generate().

   {
 
     bool success = false;
     while (!success) {
       double neutrinoEnergy = GetNeutrinoEnergy(engine);
       double neutrinoTime = GetNeutrinoTime(engine);
       success = ProcessOneNeutrino(truth, neutrinoEnergy, neutrinoTime, engine);
     }
   }

std::vector< simb::MCTruth > evgen::NueAr40CCGenerator::Generate ( CLHEP::HepRandomEngine & engine )

Definition at line 62 of file NueAr40CCGenerator.cxx.

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

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

   {
 
     std::vector<simb::MCTruth> truths;
 
     int NumberOfNu = this->GetNumberOfNeutrinos(engine);
     for (int i = 0; i < NumberOfNu; ++i) {
 
       simb::MCTruth truth;
       truth.SetOrigin(simb::kSuperNovaNeutrino);
 
       // Loop until at least one neutrino is simulated
       while (!truth.NParticles()) {
         CreateKinematicsVector(truth, engine);
       }
 
       truths.push_back(truth);
     }
 
     return truths;
   }

std::vector< double > evgen::NueAr40CCGenerator::GetIsotropicDirection ( CLHEP::HepRandomEngine & engine ) const

private

Definition at line 86 of file NueAr40CCGenerator.cxx.

Referenced by ProcessOneNeutrino().

   {
 
     CLHEP::RandFlat randFlat(engine);
 
     std::vector<double> isotropicDirection;
 
     double phi = 2 * TMath::Pi() * randFlat.fire();
     double cosTheta = 2 * randFlat.fire() - 1;
     double theta = TMath::ACos(cosTheta);
 
     // x, y, z
     isotropicDirection.push_back(cos(phi) * sin(theta));
     isotropicDirection.push_back(sin(phi) * sin(theta));
     isotropicDirection.push_back(cosTheta);
 
     return isotropicDirection;
   }

double evgen::NueAr40CCGenerator::GetNeutrinoEnergy ( CLHEP::HepRandomEngine & engine ) const

private

Definition at line 150 of file NueAr40CCGenerator.cxx.

References fEnergyProbabilityMap, fMonoenergeticNeutrinos, and fNeutrinoEnergy.

Referenced by CreateKinematicsVector().

   {
 
     if (fMonoenergeticNeutrinos) return fNeutrinoEnergy;
 
     CLHEP::RandFlat randFlat(engine);
 
     double neutrinoEnergy = 0.0;
 
     double randomNumber = randFlat.fire();
 
     // We need this to get a previous entry in the map
     std::pair<double, double> previousPair;
 
     for (std::map<double, double>::const_iterator energyProbability = fEnergyProbabilityMap.begin();
          energyProbability != fEnergyProbabilityMap.end();
          ++energyProbability) {
       if (randomNumber < energyProbability->second) {
         if (energyProbability != fEnergyProbabilityMap.begin()) {
           neutrinoEnergy =
             energyProbability->first - (energyProbability->second - randomNumber) *
                                          (energyProbability->first - previousPair.first) /
                                          (energyProbability->second - previousPair.second);
           break;
         }
         else {
           neutrinoEnergy = energyProbability->first;
           break;
         }
       }
       previousPair = *energyProbability;
     }
 
     return neutrinoEnergy;
   }

double evgen::NueAr40CCGenerator::GetNeutrinoTime ( CLHEP::HepRandomEngine & engine ) const

private

Definition at line 139 of file NueAr40CCGenerator.cxx.

References fNeutrinoTimeBegin, and fNeutrinoTimeEnd.

Referenced by CreateKinematicsVector().

   {
 
     CLHEP::RandFlat randFlat(engine);
 
     return randFlat.fire(fNeutrinoTimeBegin, fNeutrinoTimeEnd);
   }

int evgen::NueAr40CCGenerator::GetNumberOfNeutrinos ( CLHEP::HepRandomEngine & engine ) const

private

Definition at line 124 of file NueAr40CCGenerator.cxx.

References fAllowZeroNeutrinos, fNumberOfNeutrinos, and fUsePoissonDistribution.

Referenced by Generate().

   {
 
     if (fUsePoissonDistribution) {
       CLHEP::RandPoisson randPoisson(engine);
       int N = randPoisson.fire(fNumberOfNeutrinos);
       if (N == 0 && !fAllowZeroNeutrinos) N = 1;
       return N;
     }
 
     return fNumberOfNeutrinos;
   }

std::vector< double > evgen::NueAr40CCGenerator::GetUniformPosition ( CLHEP::HepRandomEngine & engine ) const

private

Definition at line 108 of file NueAr40CCGenerator.cxx.

References fActiveVolume.

Referenced by ProcessOneNeutrino().

   {
 
     CLHEP::RandFlat randFlat(engine);
 
     std::vector<double> position;
 
     position.push_back(randFlat.fire(fActiveVolume.at(0).at(0), fActiveVolume.at(1).at(0)));
     position.push_back(randFlat.fire(fActiveVolume.at(0).at(1), fActiveVolume.at(1).at(1)));
     position.push_back(randFlat.fire(fActiveVolume.at(0).at(2), fActiveVolume.at(1).at(2)));
 
     return position;
   }

void evgen::NueAr40CCGenerator::InitializeVectors ( )

private

Definition at line 223 of file NueAr40CCGenerator.cxx.

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

Referenced by NueAr40CCGenerator().

   {
 
     // The level data -- there's probably a better way to initialize
     fBranchingRatios.at(0).push_back(1.00);
     fDecayTo.at(0).push_back(1);
 
     fBranchingRatios.at(1).push_back(1.00);
     fDecayTo.at(1).push_back(4);
 
     fBranchingRatios.at(2).push_back(0.133891);
     fBranchingRatios.at(2).push_back(0.41841);
     fBranchingRatios.at(2).push_back(0.351464);
     fBranchingRatios.at(2).push_back(0.0962343);
     fDecayTo.at(2).push_back(38);
     fDecayTo.at(2).push_back(43);
     fDecayTo.at(2).push_back(59);
     fDecayTo.at(2).push_back(72);
 
     fBranchingRatios.at(3).push_back(0.391705);
     fBranchingRatios.at(3).push_back(0.460829);
     fBranchingRatios.at(3).push_back(0.147465);
     fDecayTo.at(3).push_back(61);
     fDecayTo.at(3).push_back(65);
     fDecayTo.at(3).push_back(71);
 
     fBranchingRatios.at(4).push_back(0.358974);
     fBranchingRatios.at(4).push_back(0.641026);
     fDecayTo.at(4).push_back(13);
     fDecayTo.at(4).push_back(58);
 
     fBranchingRatios.at(5).push_back(0.247808);
     fBranchingRatios.at(5).push_back(0.0468929);
     fBranchingRatios.at(5).push_back(0.213496);
     fBranchingRatios.at(5).push_back(0.11056);
     fBranchingRatios.at(5).push_back(0.381243);
     fDecayTo.at(5).push_back(40);
     fDecayTo.at(5).push_back(52);
     fDecayTo.at(5).push_back(67);
     fDecayTo.at(5).push_back(71);
     fDecayTo.at(5).push_back(72);
 
     fBranchingRatios.at(6).push_back(0.361025);
     fBranchingRatios.at(6).push_back(0.056677);
     fBranchingRatios.at(6).push_back(0.194099);
     fBranchingRatios.at(6).push_back(0.388199);
     fDecayTo.at(6).push_back(52);
     fDecayTo.at(6).push_back(55);
     fDecayTo.at(6).push_back(63);
     fDecayTo.at(6).push_back(68);
 
     fBranchingRatios.at(7).push_back(0.0300963);
     fBranchingRatios.at(7).push_back(0.0613965);
     fBranchingRatios.at(7).push_back(0.0152488);
     fBranchingRatios.at(7).push_back(0.0321027);
     fBranchingRatios.at(7).push_back(0.0513644);
     fBranchingRatios.at(7).push_back(0.0682183);
     fBranchingRatios.at(7).push_back(0.073435);
     fBranchingRatios.at(7).push_back(0.0100321);
     fBranchingRatios.at(7).push_back(0.0120385);
     fBranchingRatios.at(7).push_back(0.0922953);
     fBranchingRatios.at(7).push_back(0.152488);
     fBranchingRatios.at(7).push_back(0.401284);
     fDecayTo.at(7).push_back(17);
     fDecayTo.at(7).push_back(25);
     fDecayTo.at(7).push_back(27);
     fDecayTo.at(7).push_back(32);
     fDecayTo.at(7).push_back(49);
     fDecayTo.at(7).push_back(54);
     fDecayTo.at(7).push_back(56);
     fDecayTo.at(7).push_back(62);
     fDecayTo.at(7).push_back(63);
     fDecayTo.at(7).push_back(67);
     fDecayTo.at(7).push_back(68);
     fDecayTo.at(7).push_back(70);
 
     fBranchingRatios.at(8).push_back(0.0089877);
     fBranchingRatios.at(8).push_back(0.06386);
     fBranchingRatios.at(8).push_back(0.13245);
     fBranchingRatios.at(8).push_back(0.473037);
     fBranchingRatios.at(8).push_back(0.321665);
     fDecayTo.at(8).push_back(43);
     fDecayTo.at(8).push_back(56);
     fDecayTo.at(8).push_back(67);
     fDecayTo.at(8).push_back(70);
     fDecayTo.at(8).push_back(71);
 
     fBranchingRatios.at(9).push_back(0.16129);
     fBranchingRatios.at(9).push_back(0.193548);
     fBranchingRatios.at(9).push_back(0.645161);
     fDecayTo.at(9).push_back(37);
     fDecayTo.at(9).push_back(65);
     fDecayTo.at(9).push_back(72);
 
     fBranchingRatios.at(10).push_back(0.245826);
     fBranchingRatios.at(10).push_back(0.0816327);
     fBranchingRatios.at(10).push_back(0.20872);
     fBranchingRatios.at(10).push_back(0.463822);
     fDecayTo.at(10).push_back(40);
     fDecayTo.at(10).push_back(56);
     fDecayTo.at(10).push_back(60);
     fDecayTo.at(10).push_back(71);
 
     fBranchingRatios.at(11).push_back(0.170984);
     fBranchingRatios.at(11).push_back(0.310881);
     fBranchingRatios.at(11).push_back(0.518135);
     fDecayTo.at(11).push_back(29);
     fDecayTo.at(11).push_back(54);
     fDecayTo.at(11).push_back(66);
 
     fBranchingRatios.at(12).push_back(0.242424);
     fBranchingRatios.at(12).push_back(0.757576);
     fDecayTo.at(12).push_back(54);
     fDecayTo.at(12).push_back(62);
 
     fBranchingRatios.at(13).push_back(0.159664);
     fBranchingRatios.at(13).push_back(0.840336);
     fDecayTo.at(13).push_back(48);
     fDecayTo.at(13).push_back(58);
 
     fBranchingRatios.at(14).push_back(0.288136);
     fBranchingRatios.at(14).push_back(0.145763);
     fBranchingRatios.at(14).push_back(0.118644);
     fBranchingRatios.at(14).push_back(0.108475);
     fBranchingRatios.at(14).push_back(0.338983);
     fDecayTo.at(14).push_back(56);
     fDecayTo.at(14).push_back(66);
     fDecayTo.at(14).push_back(70);
     fDecayTo.at(14).push_back(71);
     fDecayTo.at(14).push_back(72);
 
     fBranchingRatios.at(15).push_back(0.00869263);
     fBranchingRatios.at(15).push_back(0.087274);
     fBranchingRatios.at(15).push_back(0.0973574);
     fBranchingRatios.at(15).push_back(0.288595);
     fBranchingRatios.at(15).push_back(0.347705);
     fBranchingRatios.at(15).push_back(0.170376);
     fDecayTo.at(15).push_back(40);
     fDecayTo.at(15).push_back(64);
     fDecayTo.at(15).push_back(65);
     fDecayTo.at(15).push_back(68);
     fDecayTo.at(15).push_back(70);
     fDecayTo.at(15).push_back(71);
 
     fBranchingRatios.at(16).push_back(1);
     fDecayTo.at(16).push_back(65);
 
     fBranchingRatios.at(17).push_back(0.341763);
     fBranchingRatios.at(17).push_back(0.0567327);
     fBranchingRatios.at(17).push_back(0.164046);
     fBranchingRatios.at(17).push_back(0.239234);
     fBranchingRatios.at(17).push_back(0.198223);
     fDecayTo.at(17).push_back(35);
     fDecayTo.at(17).push_back(39);
     fDecayTo.at(17).push_back(51);
     fDecayTo.at(17).push_back(67);
     fDecayTo.at(17).push_back(69);
 
     fBranchingRatios.at(18).push_back(0.106406);
     fBranchingRatios.at(18).push_back(0.254103);
     fBranchingRatios.at(18).push_back(0.0465855);
     fBranchingRatios.at(18).push_back(0.529381);
     fBranchingRatios.at(18).push_back(0.0635257);
     fDecayTo.at(18).push_back(60);
     fDecayTo.at(18).push_back(61);
     fDecayTo.at(18).push_back(63);
     fDecayTo.at(18).push_back(70);
     fDecayTo.at(18).push_back(72);
 
     fBranchingRatios.at(19).push_back(0.0905797);
     fBranchingRatios.at(19).push_back(0.228261);
     fBranchingRatios.at(19).push_back(0.181159);
     fBranchingRatios.at(19).push_back(0.137681);
     fBranchingRatios.at(19).push_back(0.362319);
     fDecayTo.at(19).push_back(43);
     fDecayTo.at(19).push_back(45);
     fDecayTo.at(19).push_back(52);
     fDecayTo.at(19).push_back(70);
     fDecayTo.at(19).push_back(71);
 
     fBranchingRatios.at(20).push_back(0.0316414);
     fBranchingRatios.at(20).push_back(0.0415293);
     fBranchingRatios.at(20).push_back(0.0362558);
     fBranchingRatios.at(20).push_back(0.0481213);
     fBranchingRatios.at(20).push_back(0.0903098);
     fBranchingRatios.at(20).push_back(0.0929466);
     fBranchingRatios.at(20).push_back(0.659196);
     fDecayTo.at(20).push_back(30);
     fDecayTo.at(20).push_back(32);
     fDecayTo.at(20).push_back(46);
     fDecayTo.at(20).push_back(61);
     fDecayTo.at(20).push_back(64);
     fDecayTo.at(20).push_back(66);
     fDecayTo.at(20).push_back(70);
 
     fBranchingRatios.at(21).push_back(0.310757);
     fBranchingRatios.at(21).push_back(0.398406);
     fBranchingRatios.at(21).push_back(0.290837);
     fDecayTo.at(21).push_back(64);
     fDecayTo.at(21).push_back(66);
     fDecayTo.at(21).push_back(70);
 
     fBranchingRatios.at(22).push_back(0.152542);
     fBranchingRatios.at(22).push_back(0.847458);
     fDecayTo.at(22).push_back(67);
     fDecayTo.at(22).push_back(71);
 
     fBranchingRatios.at(23).push_back(0.168367);
     fBranchingRatios.at(23).push_back(0.321429);
     fBranchingRatios.at(23).push_back(0.510204);
     fDecayTo.at(23).push_back(52);
     fDecayTo.at(23).push_back(70);
     fDecayTo.at(23).push_back(71);
 
     fBranchingRatios.at(24).push_back(0.0276437);
     fBranchingRatios.at(24).push_back(0.078982);
     fBranchingRatios.at(24).push_back(0.0245722);
     fBranchingRatios.at(24).push_back(0.162352);
     fBranchingRatios.at(24).push_back(0.184291);
     fBranchingRatios.at(24).push_back(0.438789);
     fBranchingRatios.at(24).push_back(0.0833699);
     fDecayTo.at(24).push_back(36);
     fDecayTo.at(24).push_back(53);
     fDecayTo.at(24).push_back(62);
     fDecayTo.at(24).push_back(64);
     fDecayTo.at(24).push_back(70);
     fDecayTo.at(24).push_back(71);
     fDecayTo.at(24).push_back(72);
 
     fBranchingRatios.at(25).push_back(0.0163934);
     fBranchingRatios.at(25).push_back(0.336276);
     fBranchingRatios.at(25).push_back(0.226986);
     fBranchingRatios.at(25).push_back(0.420345);
     fDecayTo.at(25).push_back(43);
     fDecayTo.at(25).push_back(67);
     fDecayTo.at(25).push_back(68);
     fDecayTo.at(25).push_back(70);
 
     fBranchingRatios.at(26).push_back(0.184332);
     fBranchingRatios.at(26).push_back(0.0460829);
     fBranchingRatios.at(26).push_back(0.460829);
     fBranchingRatios.at(26).push_back(0.078341);
     fBranchingRatios.at(26).push_back(0.230415);
     fDecayTo.at(26).push_back(53);
     fDecayTo.at(26).push_back(54);
     fDecayTo.at(26).push_back(60);
     fDecayTo.at(26).push_back(61);
     fDecayTo.at(26).push_back(71);
 
     fBranchingRatios.at(27).push_back(0.0147145);
     fBranchingRatios.at(27).push_back(0.0170689);
     fBranchingRatios.at(27).push_back(0.0500294);
     fBranchingRatios.at(27).push_back(0.329606);
     fBranchingRatios.at(27).push_back(0.588582);
     fDecayTo.at(27).push_back(36);
     fDecayTo.at(27).push_back(46);
     fDecayTo.at(27).push_back(56);
     fDecayTo.at(27).push_back(67);
     fDecayTo.at(27).push_back(68);
 
     fBranchingRatios.at(28).push_back(1);
     fDecayTo.at(28).push_back(70);
 
     fBranchingRatios.at(29).push_back(0.280702);
     fBranchingRatios.at(29).push_back(0.0935673);
     fBranchingRatios.at(29).push_back(0.0409357);
     fBranchingRatios.at(29).push_back(0.584795);
     fDecayTo.at(29).push_back(63);
     fDecayTo.at(29).push_back(66);
     fDecayTo.at(29).push_back(68);
     fDecayTo.at(29).push_back(70);
 
     fBranchingRatios.at(30).push_back(0.393701);
     fBranchingRatios.at(30).push_back(0.283465);
     fBranchingRatios.at(30).push_back(0.188976);
     fBranchingRatios.at(30).push_back(0.133858);
     fDecayTo.at(30).push_back(61);
     fDecayTo.at(30).push_back(67);
     fDecayTo.at(30).push_back(71);
     fDecayTo.at(30).push_back(72);
 
     fBranchingRatios.at(31).push_back(0.0477737);
     fBranchingRatios.at(31).push_back(0.194805);
     fBranchingRatios.at(31).push_back(0.0245826);
     fBranchingRatios.at(31).push_back(0.269017);
     fBranchingRatios.at(31).push_back(0.463822);
     fDecayTo.at(31).push_back(45);
     fDecayTo.at(31).push_back(60);
     fDecayTo.at(31).push_back(65);
     fDecayTo.at(31).push_back(71);
     fDecayTo.at(31).push_back(72);
 
     fBranchingRatios.at(32).push_back(0.105769);
     fBranchingRatios.at(32).push_back(0.129808);
     fBranchingRatios.at(32).push_back(0.0528846);
     fBranchingRatios.at(32).push_back(0.480769);
     fBranchingRatios.at(32).push_back(0.230769);
     fDecayTo.at(32).push_back(46);
     fDecayTo.at(32).push_back(56);
     fDecayTo.at(32).push_back(66);
     fDecayTo.at(32).push_back(70);
     fDecayTo.at(32).push_back(71);
 
     fBranchingRatios.at(33).push_back(0.21875);
     fBranchingRatios.at(33).push_back(0.78125);
     fDecayTo.at(33).push_back(67);
     fDecayTo.at(33).push_back(71);
 
     fBranchingRatios.at(34).push_back(0.102011);
     fBranchingRatios.at(34).push_back(0.179598);
     fBranchingRatios.at(34).push_back(0.718391);
     fDecayTo.at(34).push_back(43);
     fDecayTo.at(34).push_back(61);
     fDecayTo.at(34).push_back(66);
 
     fBranchingRatios.at(35).push_back(0.00826446);
     fBranchingRatios.at(35).push_back(0.393546);
     fBranchingRatios.at(35).push_back(0.334514);
     fBranchingRatios.at(35).push_back(0.263676);
     fDecayTo.at(35).push_back(64);
     fDecayTo.at(35).push_back(67);
     fDecayTo.at(35).push_back(68);
     fDecayTo.at(35).push_back(70);
 
     fBranchingRatios.at(36).push_back(0.056338);
     fBranchingRatios.at(36).push_back(0.704225);
     fBranchingRatios.at(36).push_back(0.239437);
     fDecayTo.at(36).push_back(51);
     fDecayTo.at(36).push_back(70);
     fDecayTo.at(36).push_back(71);
 
     fBranchingRatios.at(37).push_back(0.21875);
     fBranchingRatios.at(37).push_back(0.78125);
     fDecayTo.at(37).push_back(67);
     fDecayTo.at(37).push_back(70);
 
     fBranchingRatios.at(38).push_back(0.181818);
     fBranchingRatios.at(38).push_back(0.757576);
     fBranchingRatios.at(38).push_back(0.0606061);
     fDecayTo.at(38).push_back(66);
     fDecayTo.at(38).push_back(71);
     fDecayTo.at(38).push_back(72);
 
     fBranchingRatios.at(39).push_back(0.157258);
     fBranchingRatios.at(39).push_back(0.403226);
     fBranchingRatios.at(39).push_back(0.237903);
     fBranchingRatios.at(39).push_back(0.201613);
     fDecayTo.at(39).push_back(62);
     fDecayTo.at(39).push_back(70);
     fDecayTo.at(39).push_back(71);
     fDecayTo.at(39).push_back(72);
 
     fBranchingRatios.at(40).push_back(0.0740741);
     fBranchingRatios.at(40).push_back(0.925926);
     fDecayTo.at(40).push_back(52);
     fDecayTo.at(40).push_back(72);
 
     fBranchingRatios.at(41).push_back(0.0535714);
     fBranchingRatios.at(41).push_back(0.35119);
     fBranchingRatios.at(41).push_back(0.595238);
     fDecayTo.at(41).push_back(67);
     fDecayTo.at(41).push_back(68);
     fDecayTo.at(41).push_back(70);
 
     fBranchingRatios.at(42).push_back(0.00816803);
     fBranchingRatios.at(42).push_back(0.0583431);
     fBranchingRatios.at(42).push_back(0.350058);
     fBranchingRatios.at(42).push_back(0.583431);
     fDecayTo.at(42).push_back(49);
     fDecayTo.at(42).push_back(62);
     fDecayTo.at(42).push_back(71);
     fDecayTo.at(42).push_back(72);
 
     fBranchingRatios.at(43).push_back(0.0961538);
     fBranchingRatios.at(43).push_back(0.423077);
     fBranchingRatios.at(43).push_back(0.480769);
     fDecayTo.at(43).push_back(66);
     fDecayTo.at(43).push_back(67);
     fDecayTo.at(43).push_back(68);
 
     fBranchingRatios.at(44).push_back(0.450549);
     fBranchingRatios.at(44).push_back(0.549451);
     fDecayTo.at(44).push_back(69);
     fDecayTo.at(44).push_back(72);
 
     fBranchingRatios.at(45).push_back(0.469925);
     fBranchingRatios.at(45).push_back(0.0836466);
     fBranchingRatios.at(45).push_back(0.446429);
     fDecayTo.at(45).push_back(61);
     fDecayTo.at(45).push_back(65);
     fDecayTo.at(45).push_back(72);
 
     fBranchingRatios.at(46).push_back(0.0408163);
     fBranchingRatios.at(46).push_back(0.510204);
     fBranchingRatios.at(46).push_back(0.44898);
     fDecayTo.at(46).push_back(67);
     fDecayTo.at(46).push_back(70);
     fDecayTo.at(46).push_back(71);
 
     fBranchingRatios.at(47).push_back(1);
     fDecayTo.at(47).push_back(72);
 
     fBranchingRatios.at(48).push_back(0.641026);
     fBranchingRatios.at(48).push_back(0.358974);
     fDecayTo.at(48).push_back(58);
     fDecayTo.at(48).push_back(69);
 
     fBranchingRatios.at(49).push_back(0.0458015);
     fBranchingRatios.at(49).push_back(0.954198);
     fDecayTo.at(49).push_back(66);
     fDecayTo.at(49).push_back(70);
 
     fBranchingRatios.at(50).push_back(0.401639);
     fBranchingRatios.at(50).push_back(0.188525);
     fBranchingRatios.at(50).push_back(0.409836);
     fDecayTo.at(50).push_back(61);
     fDecayTo.at(50).push_back(69);
     fDecayTo.at(50).push_back(72);
 
     fBranchingRatios.at(51).push_back(0.0188679);
     fBranchingRatios.at(51).push_back(0.173585);
     fBranchingRatios.at(51).push_back(0.754717);
     fBranchingRatios.at(51).push_back(0.0528302);
     fDecayTo.at(51).push_back(61);
     fDecayTo.at(51).push_back(67);
     fDecayTo.at(51).push_back(71);
     fDecayTo.at(51).push_back(72);
 
     fBranchingRatios.at(52).push_back(0.0100849);
     fBranchingRatios.at(52).push_back(0.00796178);
     fBranchingRatios.at(52).push_back(0.530786);
     fBranchingRatios.at(52).push_back(0.451168);
     fDecayTo.at(52).push_back(59);
     fDecayTo.at(52).push_back(68);
     fDecayTo.at(52).push_back(70);
     fDecayTo.at(52).push_back(71);
 
     fBranchingRatios.at(53).push_back(0.0379902);
     fBranchingRatios.at(53).push_back(0.0490196);
     fBranchingRatios.at(53).push_back(0.612745);
     fBranchingRatios.at(53).push_back(0.300245);
     fDecayTo.at(53).push_back(67);
     fDecayTo.at(53).push_back(70);
     fDecayTo.at(53).push_back(71);
     fDecayTo.at(53).push_back(72);
 
     fBranchingRatios.at(54).push_back(0.106557);
     fBranchingRatios.at(54).push_back(0.819672);
     fBranchingRatios.at(54).push_back(0.0737705);
     fDecayTo.at(54).push_back(59);
     fDecayTo.at(54).push_back(68);
     fDecayTo.at(54).push_back(70);
 
     fBranchingRatios.at(55).push_back(0.699301);
     fBranchingRatios.at(55).push_back(0.300699);
     fDecayTo.at(55).push_back(64);
     fDecayTo.at(55).push_back(70);
 
     fBranchingRatios.at(56).push_back(1);
     fDecayTo.at(56).push_back(71);
 
     fBranchingRatios.at(57).push_back(1);
     fDecayTo.at(57).push_back(72);
 
     fBranchingRatios.at(58).push_back(0.888099);
     fBranchingRatios.at(58).push_back(0.111901);
     fDecayTo.at(58).push_back(69);
     fDecayTo.at(58).push_back(72);
 
     fBranchingRatios.at(59).push_back(0.00647298);
     fBranchingRatios.at(59).push_back(0.752672);
     fBranchingRatios.at(59).push_back(0.165588);
     fBranchingRatios.at(59).push_back(0.0752672);
     fDecayTo.at(59).push_back(65);
     fDecayTo.at(59).push_back(70);
     fDecayTo.at(59).push_back(71);
     fDecayTo.at(59).push_back(72);
 
     fBranchingRatios.at(60).push_back(0.0708556);
     fBranchingRatios.at(60).push_back(0.668449);
     fBranchingRatios.at(60).push_back(0.260695);
     fDecayTo.at(60).push_back(65);
     fDecayTo.at(60).push_back(71);
     fDecayTo.at(60).push_back(72);
 
     fBranchingRatios.at(61).push_back(0.166667);
     fBranchingRatios.at(61).push_back(0.833333);
     fDecayTo.at(61).push_back(69);
     fDecayTo.at(61).push_back(72);
 
     fBranchingRatios.at(62).push_back(0.0898551);
     fBranchingRatios.at(62).push_back(0.57971);
     fBranchingRatios.at(62).push_back(0.330435);
     fDecayTo.at(62).push_back(67);
     fDecayTo.at(62).push_back(68);
     fDecayTo.at(62).push_back(70);
 
     fBranchingRatios.at(63).push_back(0.813008);
     fBranchingRatios.at(63).push_back(0.186992);
     fDecayTo.at(63).push_back(71);
     fDecayTo.at(63).push_back(72);
 
     fBranchingRatios.at(64).push_back(0.29078);
     fBranchingRatios.at(64).push_back(0.70922);
     fDecayTo.at(64).push_back(70);
     fDecayTo.at(64).push_back(71);
 
     fBranchingRatios.at(65).push_back(0.05);
     fBranchingRatios.at(65).push_back(0.08);
     fBranchingRatios.at(65).push_back(0.5);
     fBranchingRatios.at(65).push_back(0.37);
     fDecayTo.at(65).push_back(69);
     fDecayTo.at(65).push_back(70);
     fDecayTo.at(65).push_back(71);
     fDecayTo.at(65).push_back(72);
 
     fBranchingRatios.at(66).push_back(0.398406);
     fBranchingRatios.at(66).push_back(0.310757);
     fBranchingRatios.at(66).push_back(0.290837);
     fDecayTo.at(66).push_back(70);
     fDecayTo.at(66).push_back(71);
     fDecayTo.at(66).push_back(72);
 
     fBranchingRatios.at(67).push_back(0.819672);
     fBranchingRatios.at(67).push_back(0.180328);
     fDecayTo.at(67).push_back(70);
     fDecayTo.at(67).push_back(71);
 
     fBranchingRatios.at(68).push_back(0.186992);
     fBranchingRatios.at(68).push_back(0.813008);
     fDecayTo.at(68).push_back(70);
     fDecayTo.at(68).push_back(71);
 
     fBranchingRatios.at(69).push_back(1);
     fDecayTo.at(69).push_back(72);
 
     fBranchingRatios.at(70).push_back(1);
     fDecayTo.at(70).push_back(71);
 
     fBranchingRatios.at(71).push_back(1);
     fDecayTo.at(71).push_back(72);
 
     // Info from table VII (http://prc.aps.org/pdf/PRC/v58/i6/p3677_1)
     // in MeV
     double startEnergyLevels[] = {2.281, 2.752, 2.937, 3.143, 3.334, 3.569, 3.652,
                                   3.786, 3.861, 4.067, 4.111, 4.267, 4.364, 4.522,
                                   4.655, 4.825, 5.017, 5.080, 5.223, 5.696, 6.006};
 
     fStartEnergyLevels.insert(
       fStartEnergyLevels.end(), &startEnergyLevels[0], &startEnergyLevels[fNumberOfStartLevels]);
 
     double b[] = {0.9,  1.5,  0.11, 0.06, 0.04, 0.01, 0.16, 0.26, 0.01, 0.05, 0.11,
                   0.29, 3.84, 0.31, 0.38, 0.47, 0.36, 0.23, 0.03, 0.11, 0.13};
 
     fB.insert(fB.end(), &b[0], &b[fNumberOfStartLevels]);
 
     double energyLevels[] = {
       7.4724,   6.2270,   5.06347,  4.99294,  4.8756,   4.87255,  4.78865,  4.744093,  4.53706,
       4.47299,  4.41936,  4.39588,  4.3840,   4.3656,   4.28052,  4.25362,  4.21307,   4.18003,
       4.14901,  4.11084,  4.10446,  4.02035,  3.92390,  3.88792,  3.86866,  3.840228,  3.82143,
       3.79757,  3.76779,  3.73848,  3.663739, 3.62995,  3.59924,  3.55697,  3.48621,   3.439144,
       3.41434,  3.39363,  3.36803,  3.22867,  3.15381,  3.14644,  3.12836,  3.109721,  3.1002,
       3.02795,  2.98587,  2.9508,   2.87901,  2.80788,  2.7874,   2.786644, 2.75672,   2.74691,
       2.730372, 2.625990, 2.57593,  2.558,    2.54277,  2.419171, 2.397165, 2.290493,  2.289871,
       2.26040,  2.103668, 2.069809, 2.047354, 1.959068, 1.643639, 0.891398, 0.8001427, 0.0298299,
       0.0};
 
     fEnergyLevels.insert(fEnergyLevels.end(), &energyLevels[0], &energyLevels[fNumberOfLevels]);
 
     // I have to put checks that the arrays really
     // have as many elements as they should
   }

bool evgen::NueAr40CCGenerator::ProcessOneNeutrino	(	simb::MCTruth &	truth,
		double	neutrinoEnergy,
		double	neutrinoTime,
		CLHEP::HepRandomEngine &	engine
	)		const

private

Definition at line 813 of file NueAr40CCGenerator.cxx.

References util::abs(), 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().

   {
 
     CLHEP::RandFlat randFlat(engine);
 
     int highestLevel = 0;
     std::vector<double> levelCrossSections = CalculateCrossSections(neutrinoEnergy, highestLevel);
 
     double totalCrossSection = 0;
     // Calculating total cross section
     for (std::vector<double>::iterator crossSection = levelCrossSections.begin();
          crossSection != levelCrossSections.end();
          ++crossSection)
       totalCrossSection += *crossSection;
 
     if (totalCrossSection == 0) return false;
 
     std::vector<double> startLevelProbabilities;
     // Calculating each level's probability
     for (std::vector<double>::iterator crossSection = levelCrossSections.begin();
          crossSection != levelCrossSections.end();
          ++crossSection)
       startLevelProbabilities.push_back((*crossSection) / totalCrossSection);
 
     double randomNumber = randFlat.fire();
     double tprob = 0;
     int chosenStartLevel = -1;
     // Picking a starting level
     for (int level = 0; level < highestLevel; ++level) {
       if (randomNumber < (startLevelProbabilities.at(level) + tprob)) {
         chosenStartLevel = level;
         break;
       }
       tprob += startLevelProbabilities.at(level);
     }
 
     int lastLevel = -1;
     int level = -1;
 
     // Time delays
     // Seems like this is not implemented yet
     std::vector<double> levelDelay;
     for (int n = 0; n < fNumberOfLevels; ++n)
       levelDelay.push_back(0.0);
 
     // The highest n for which fStartEnergyLevels.at(chosenStartLevel)
     // is higher than fEnergyLevels.at(n)
     int highestHigher = 0;
     // The lowest n for which fStartEnergyLevels.at(chosenStartLevel)
     // is lower than fEnergyLevels.at(n)
     int lowestLower = 0;
 
     // Finding lowestLower and highestHigher
     for (int n = 0; n < fNumberOfLevels; ++n) {
       if (fStartEnergyLevels.at(chosenStartLevel) < fEnergyLevels.at(n)) lowestLower = n;
       if (fStartEnergyLevels.at(chosenStartLevel) > fEnergyLevels.at(n)) {
         highestHigher = n;
         break;
       }
     }
 
     if (std::abs(fStartEnergyLevels.at(chosenStartLevel) - fEnergyLevels.at(lowestLower)) <
         std::abs(fStartEnergyLevels.at(chosenStartLevel) - fEnergyLevels.at(highestHigher))) {
       lastLevel = lowestLower;
       level = lowestLower;
     }
     // If the chosen start level energy is closest to the lowest level energy
     // that it's lower than than the highest level energy
     // that it's higher than, it starts at the level of
     // the lowest level energy that it's lower than
 
     if (std::abs(fStartEnergyLevels.at(chosenStartLevel) - fEnergyLevels.at(highestHigher)) <
         std::abs(fStartEnergyLevels.at(chosenStartLevel) - fEnergyLevels.at(lowestLower))) {
       lastLevel = highestHigher;
       level = highestHigher;
     }
     // If the chosen start level energy is closest to the highest level energy
     // that it's higher than than the lowest level energy that it's lower than,
     // it starts at the level of the highest level energy that it's higher than
 
     std::vector<double> vertex = GetUniformPosition(engine);
 
     std::vector<double> neutrinoDirection = GetIsotropicDirection(engine);
 
     //
     // Add the first particle (the neutrino) to the truth list...
     //
 
     // NOTE: all of the values below calculated for the neutrino should be checked!
 
     // Primary particles have negative IDs
     int neutrinoTrackID = -1 * (truth.NParticles() + 1);
     std::string primary("primary");
 
     int nuePDG = 12;
     double neutrinoP = neutrinoEnergy / 1000.0; // use GeV...
     double neutrinoPx = neutrinoDirection.at(0) * neutrinoP;
     double neutrinoPy = neutrinoDirection.at(1) * neutrinoP;
     double neutrinoPz = neutrinoDirection.at(2) * neutrinoP;
 
     // Create the neutrino:
     //   * set the mother to -1 since it is primary
     //   * set the mass to something < 0 so that the constructor looks up the mass from the pdg tables
     //   * set the status bit to 0 so that geant doesn't waste any CPU tracking the neutrino
     simb::MCParticle neutrino(neutrinoTrackID, nuePDG, primary, -1, -1, 0);
 
     TLorentzVector neutrinoPosition(vertex.at(0), vertex.at(1), vertex.at(2), neutrinoTime);
     TLorentzVector neutrinoMomentum(neutrinoPx, neutrinoPy, neutrinoPz, neutrinoEnergy / 1000.0);
     neutrino.AddTrajectoryPoint(neutrinoPosition, neutrinoMomentum);
 
     truth.Add(neutrino);
 
     //
     // Adding the electron to truth
     //
 
     // In MeV
     double electronEnergy = neutrinoEnergy - (fEnergyLevels.at(lastLevel) + 1.5);
     double electronEnergyGeV = electronEnergy / 1000.0;
 
     double electronM = 0.000511;
     double electronP = std::sqrt(std::pow(electronEnergyGeV, 2) - std::pow(electronM, 2));
 
     // For the moment, choose an isotropic direction for the electron.
     std::vector<double> electronDirection = GetIsotropicDirection(engine);
     double electronPx = electronDirection.at(0) * electronP;
     double electronPy = electronDirection.at(1) * electronP;
     double electronPz = electronDirection.at(2) * electronP;
 
     // Primary particles have negative IDs
     int trackID = -1 * (truth.NParticles() + 1);
     int electronPDG = 11;
     simb::MCParticle electron(trackID, electronPDG, primary);
 
     TLorentzVector electronPosition(vertex.at(0), vertex.at(1), vertex.at(2), neutrinoTime);
     TLorentzVector electronMomentum(electronPx, electronPy, electronPz, electronEnergyGeV);
     electron.AddTrajectoryPoint(electronPosition, electronMomentum);
 
     truth.Add(electron);
 
     double ttime = neutrinoTime;
     int noMoreDecay = 0;
 
     // Level loop
     int groundLevel = fNumberOfLevels - 1;
     while (level != groundLevel) {
 
       double rl = randFlat.fire();
 
       int decayNum = 0;
 
       tprob = 0; // Used this variable above for cross section
 
       // Decay loop
       for (unsigned int iLevel = 0; iLevel < fBranchingRatios.at(level).size(); ++iLevel) {
 
         if (rl < (fBranchingRatios.at(level).at(iLevel) + tprob)) {
 
           // We have a decay
 
           level = fDecayTo.at(level).at(decayNum);
 
           // Really should be using a map
           double gammaEnergy = fEnergyLevels.at(lastLevel) - fEnergyLevels.at(level);
           double gammaEnergyGeV = gammaEnergy / 1000;
 
           std::vector<double> gammaDirection = GetIsotropicDirection(engine);
           double gammaPx = gammaDirection.at(0) * gammaEnergyGeV;
           double gammaPy = gammaDirection.at(1) * gammaEnergyGeV;
           double gammaPz = gammaDirection.at(2) * gammaEnergyGeV;
 
           //double gammaM  = 0.0; // unused
 
           double gammaTime =
             (-TMath::Log(randFlat.fire()) / (1 / (levelDelay.at(lastLevel)))) + ttime;
 
           // Adding the gamma to truth
 
           // Primary particles have negative IDs
           trackID = -1 * (truth.NParticles() + 1);
           int gammaPDG = 22;
           simb::MCParticle gamma(trackID, gammaPDG, primary); // NOTE: should the gammas be primary?
 
           TLorentzVector gammaPosition(vertex.at(0), vertex.at(1), vertex.at(2), neutrinoTime);
           TLorentzVector gammaMomentum(gammaPx, gammaPy, gammaPz, gammaEnergyGeV);
           gamma.AddTrajectoryPoint(gammaPosition, gammaMomentum);
 
           truth.Add(gamma);
 
           lastLevel = level;
           ttime = gammaTime;
 
           break;
         }
 
         if ((tprob + fBranchingRatios.at(level).at(iLevel)) > 1) {
           std::cout << "(tprob + *iLevel) > 1" << std::endl;
           noMoreDecay = 1; // If it doesn't do any more gamma decay
           break;
         }
 
         ++decayNum;
         tprob += fBranchingRatios.at(level).at(iLevel);
 
       } // End of decay loop
 
       if (noMoreDecay == 1) break;
 
     } // End of level loop
 
     if (level != 72) {
       std::cout << "level != 72" << std::endl;
       std::cout << "level  = " << level << std::endl;
     }
 
     // Set the neutrino for the MCTruth object:
     // NOTE: currently these parameters are all pretty much a guess...
     truth.SetNeutrino(
       simb::kCC,
       simb::kQE,   // not sure what the mode should be here, assumed that these are all QE???
       simb::kCCQE, // not sure what the int_type should be either...
       0,           // target is AR40? not sure how to specify that...
       0,           // nucleon pdg ???
       0,           // quark pdg ???
       -1.0,        // W ??? - not sure how to calculate this from the above
       -1.0,        // X ??? - not sure how to calculate this from the above
       -1.0,        // Y ??? - not sure how to calculate this from the above
       -1.0);       // Qsqr ??? - not sure how to calculate this from the above
 
     return true;
   }

void evgen::NueAr40CCGenerator::ReadNeutrinoSpectrum ( )

private

Definition at line 188 of file NueAr40CCGenerator.cxx.

References fEnergyProbabilityMap, and fEnergySpectrumFileName.

Referenced by NueAr40CCGenerator().

   {
 
     cet::search_path searchPath("FW_SEARCH_PATH");
     std::string directoryName = "SupernovaNeutrinos/" + fEnergySpectrumFileName;
 
     std::string fullName;
     searchPath.find_file(directoryName, fullName);
 
     if (fullName.empty())
       throw cet::exception("NueAr40CCGenerator")
         << "Cannot find file with neutrino energy spectrum " << fullName << '\n';
 
     TFile energySpectrumFile(fullName.c_str(), "READ");
 
     std::string energySpectrumGraphName = "NueSpectrum";
     TGraph* energySpectrumGraph =
       dynamic_cast<TGraph*>(energySpectrumFile.Get(energySpectrumGraphName.c_str()));
 
     double integral = 0.0;
     int numberOfPoints = energySpectrumGraph->GetN();
     double* energyValues = energySpectrumGraph->GetX();
     double* fluxValues = energySpectrumGraph->GetY();
     for (int point = 0; point < numberOfPoints; ++point)
       integral += fluxValues[point];
 
     double probability = 0.0;
     for (int point = 0; point < numberOfPoints; ++point) {
       probability += fluxValues[point] / integral;
       fEnergyProbabilityMap.insert(std::make_pair(energyValues[point], probability));
     }
   }