NueAr40CCGenerator.h

Go to the documentation of this file.

//=============================================================================
// NueAr40CCGenerator.h 
//
// Gleb Sinev, Duke, 2015
// Supernova neutrino generator that simulates nu_e-Ar40 CC interactions
// and produces a text file with kinematics information,
// which can be used as input for LArSoft
// Based on a generator written by AJ Roeth
//=============================================================================

// C++ includes
#include <vector>
#include <map>
#include <iostream>
#include <string>

// Framework includes
#include "fhiclcpp/ParameterSet.h"

// nutools includes
#include "nusimdata/SimulationBase/MCTruth.h"

#include "CLHEP/Random/RandomEngine.h"

namespace evgen {

  class NueAr40CCGenerator {
  
    public:
  
      // Constructor
      NueAr40CCGenerator(fhicl::ParameterSet const& parameterSet);
  
      // Simulate interactions, produce plots
      std::vector<simb::MCTruth> Generate(CLHEP::HepRandomEngine& engine);
  
    private:
  
      // Generate a direction vector isotropically
      std::vector< double > GetIsotropicDirection
                                    (CLHEP::HepRandomEngine& engine) const;
  
      // Generate a position distributed uniformly inside the active volume
      std::vector< double > GetUniformPosition
                                    (CLHEP::HepRandomEngine& engine) const;
  
      // Decide how many neutrinos to generate
      int GetNumberOfNeutrinos(CLHEP::HepRandomEngine& engine) const;

      // Generate a neutrino time distributed uniformly
      // from fNeutrinoTimeBegin to fNeutrinoTimeEnd 
      double GetNeutrinoTime(CLHEP::HepRandomEngine& engine) const;
  
      // Use fEnergyProbabilityMap to sample one energy value
      // or return a constant value
      double GetNeutrinoEnergy(CLHEP::HepRandomEngine& engine) const;
  
      // Read a ROOT file with a TGraph  and fill fEnergyProbabilityMap
      // Assume that the first point in TGraph is { 0, 0 }
      void ReadNeutrinoSpectrum();
  
      // Initialize vectors with branching ratios, probabilities,
      // and other quantities required to calculate
      // number and energy of gammas produced
      void InitializeVectors();
      
      // Simulate particles
      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; 
  
      // Assume that the first entry is { 0, 0 }
      std::map< double, double > fEnergyProbabilityMap;
  
      int fNumberOfLevels;
      int fNumberOfStartLevels;
  
      // Vectors containing information about levels
      std::vector< std::vector< double > > fBranchingRatios;
      std::vector< std::vector< int    > > fDecayTo;
      std::vector< double > fStartEnergyLevels;
      std::vector< double > fB;
      std::vector< double > fEnergyLevels;
  
      // Generate monoenergetic neutrinos if this variable is set to true
      bool fMonoenergeticNeutrinos;
      // Energy of monoenergetic neutrinos
      double fNeutrinoEnergy;

      // Otherwise sample the spectrum
      std::string fEnergySpectrumFileName;
  
      // Number of neutrinos is distributed according 
      // to Poisson distribution if this is true
      bool fUsePoissonDistribution;

      // Allow zero neutrinos to be created if that is what the randonmly
      // generated number-of-nu ends up being.
      bool fAllowZeroNeutrinos;
      
      // Average or exact number of neutrinos generated
      int fNumberOfNeutrinos;
     
      // Generate neutrinos uniformly in a time interval
      double fNeutrinoTimeBegin;
      double fNeutrinoTimeEnd;
  
      // Two opposite vertices { { x0, y0, z0 }, { x1, y1, z1 } } 
      // of the active volume box, such that x0 < x1, y0 < y1, z0 < z1
      std::vector< std::vector < double > > fActiveVolume;
  
  };

}