//////////////////////////////////////////////////////////////////////// /// \class RAT::ReactorGen /// \author Nuno Barros -- contact person /// \date 13-Dec-2011 /// /// \brief Specific reactor neutrino generator with configurable sub-generators for time and position. /// /// REVISION HISTORY: /// - 28-Oct-2015 : Nuno Barros /// - Initial revision /// - 17-Oct-2016 : Nuno Barros /// - Major overhaul to make the generator the finalized, production /// ready version. /// - 07-Feb-2017 : Nuno Barros /// - Cleaned up low level generators (inherited from GLG4Gen) /// - 11-Dec-2017 : Ian Coulter /// - Move LLA and Distance tools to utilities /// - 11-Oct-2019 : Nuno Barros /// - Adapted code to use the cross section manager /// /// - 16-Jun-2020 : Charlie Mills /// - Added user option for reactor core power scaling /// //////////////////////////////////////////////////////////////////////// #ifndef __RAT_ReactorGen_ #define __RAT_ReactorGen_ /// RAT includes #include #include #include /// Geant4 includes #include #include #include /// STL #include #include class TH1D; class TH2D; using std::string; namespace RAT { /// Forward declarations inside namespace class VertexGen_IBD; // This is the C++ equivalent of a RATDB entry for a reactor class Reactor { public: struct CoreData { G4String reactor_name; G4double core_power; G4String core_spectrum; G4double core_distance; // from the detector G4ThreeVector core_LLA_coord; G4ThreeVector core_direction; bool is_active; }; Reactor(G4String name) : fName(name) { }; virtual ~Reactor() { }; void SetName(G4String name) {fName = name;}; G4String GetName() {return fName;}; void SetCore(CoreData &core, const unsigned int idx) {fCores[idx] = core;}; CoreData &GetCore(const unsigned int idx) { return fCores.at(idx);}; size_t GetNumCores() const {return fCores.size();}; private: Reactor(); std::map fCores; G4String fName; }; class ReactorSpectrum { public: enum SpectrumType{unknown = -1, polynomial=0, parameterized, shape}; virtual ~ReactorSpectrum(); ReactorSpectrum() : _e_min(0.0), _e_max(0.0), _spec_type(unknown), _spectrum(0),_spectrum_norm(0.0),_name("") {} ReactorSpectrum(G4String name,SpectrumType stype, double emin, double emax, std::vector spec_e, std::vector spec_f, std::vector isotopes, std::vector composition); void BuildIsotopeSpectrum(G4String isotope); static G4String SpectrumTypeToName(SpectrumType tp) { switch(tp) { case polynomial: return "polynomial"; break; case shape: return "shape"; break; case parameterized: return "parameterized"; break; case unknown: default: return "unknown"; break; } } static SpectrumType SpectrumNameToType(G4String nm) { if (nm == "polynomial") { return polynomial; } else if (nm == "parameterized") { return parameterized; } else if (nm == "shape") { return shape; } else { return unknown; } } double GetEmin() const {return _e_min;} void SetEmin(double emin) {_e_min = emin;} double GetEmax() const {return _e_max;} void SetEmax(double emax) {_e_max = emax;} SpectrumType GetSpecType() const {return _spec_type;} void SetSpecType(SpectrumType stype) {_spec_type = stype;} double GetSpectrumNorm() {return _spectrum_norm;} TH1D* GetSpectrum() {return _spectrum;} private: void BuildSpectrumFromIsotopes(std::vector &isotopes, std::vector &composition); TH1D * BuildSpectrumFromArrays(std::vector spec_e, std::vector spec_f); double _e_min; double _e_max; SpectrumType _spec_type; TH1D* _spectrum; double _spectrum_norm; G4String _name; }; class ReactorGen: public GLG4Gen { public: /** Default constructor*/ ReactorGen(); /** Destructor */ virtual ~ReactorGen(); /** * \brief Generates an event * * The workhorse of the generator. It calls all other methods and produces an event which is then returned by reference. * Implements GLG4Gen::GenerateEvent * * \param[out] event Generated event. */ virtual void GenerateEvent(G4Event *event); /** * \brief Resets the event time. * * Implements GLG4Gen::ResetTime * * \param[in] offset to be applied to the next event time. */ virtual void ResetTime(double offset=0.0); /** Defines repeatability of the generator. One can invoke multiple instances * of the generator to produce different sets of reactor groups. * */ virtual bool IsRepeatable() const { return true; }; /// /// State control. Controls the high level part of the methods /// /** * State setter. * * Method that sets the state for this generator instance. * The state string has the form * [position_gen]:[time_gen]:[vertex_gen]:[reactor_list] where * - [position_gen] is the position generator : default 'fill' * - [time_gen] is the time generator : default 'poisson' * - [vertex_gen] is the vertex generator : default 'ibd' * - [reactor_list] is the index out of REACTOR_LIST that should be used. Defaults 'DEFAULT' * * \param[in] state Encoded string with the generator configuration. */ virtual void SetState(G4String state); /** State getter. Returns an encoded string with the present state. */ virtual G4String GetState() const; /** * Sets the state of the internal time generator. * * \param[in] "" New state for the time generator * \sa GLG4TimeGen * */ virtual void SetTimeState(G4String state); /** Return the present state of the internal time generator. */ virtual G4String GetTimeState() const; /** * Set the state of the vertex generator. Implements GLG4Gen::SetVertexState. * \attention This is a dummy method. The vertex generator has its own messenger class to configure. States passed by here are not checked and directly forwarded to the vertex generator. * * @param "" String containing the new state. */ virtual void SetVertexState(G4String state); /** Getter for the vertex generator state. */ virtual G4String GetVertexState() const; /** Setter for state of the position generator. No checks performed here. */ virtual void SetPosState(G4String state); /** Getter for state of the position generator. */ virtual G4String GetPosState() const; /// -- This is when things should be loaded from the database. /// Still hasn't been implemented. virtual void BeginOfRun(); virtual void EndOfRun(); // // Auxiliary methods. // /** * Return the presently set Nu direction * @return vector with the present nu direction on the generator. */ G4ThreeVector GetNuDirection() { return fNuDir;}; /** * Sets the default neutrino generator direction. * @param[in] dir direction. */ void SetNuDirection(G4String dir); /** * Sets the time generator to use. * @param[in] generator Name of the generator to be used. * @param[in] force Force the instantiation of the generator (default: false). */ void SetTimeGenerator(G4String generator, bool force = false); /** * Returns the time generator to use. * @return string describing the time generator. */ G4String GetTimeGenerator() {return fTimeGenName;}; /** * Sets the position generator to use. * @param[in] generator Name of the generator to be used. * @param[in] force Force the instantiation of the generator (default: false). */ void SetPosGenerator(G4String generator, bool force = false); /** * Returns the position generator in use. * @return string describing the time generator. */ G4String GetPosGenerator() {return fPosGenName;}; /** * Returns the position generator in use. * @return string describing the time generator. */ G4String GetVertexGenerator() {return fVertexGenName;}; /** * Sets the scale for the flux being generated (in units of SSM). * @param[in] scale Number of SSMs to scale the flux (default 1.0). */ void CalculateTotalRate(); /** * Returns the scale of the flux in the generator in units of SSM. * @return the scale of the flux. */ G4double GetFluxScale() {return fFluxScale;}; void SetFluxScale(double val) {fFluxScale = val;}; /** Set the core power scaling flag */ void SetCorePowerScaling() {applyCorePowerscaling = true;}; /** * Returns the event rate for this generator. * Depends on the flux, its scale, the nu flavor and cross section. * * @return Event rate (in Hz). */ G4double GetEvtRate() {return fEvtRate;}; /** * Sets the reactor list to be loaded. */ void SetReactorList(G4String reactor_list) {fReactorListName = reactor_list;} const G4String GetReactorList() {return fReactorListName;} protected: /** * Sets the event rate in Hz. Implemented only for debug. * * @return Event rate (in Hz). */ void SetEvtRate(G4double rate) {fEvtRate = rate;}; /** * Triggers the generation of a new incoming neutrino direction. * * This method triggers the generation of a new direction for the neutrino, * considering the given options for reactors. */ void GenerateNuDirection(); static TVector3 ReactorDirection(double &latitude, double &longitude); static const G4ThreeVector SNO_LLA_coord_; static const G4ThreeVector SNO_ECEF_coord_; void LoadSpectrum(G4String name); /** \brief Auxiliary function that generates a random number for a given interval.*/ G4double GetRandomNumber(const double rmin, const double rmax); /** \brief Manages the state passed by GSim concerning the generator "flux:flavor". */ G4String fStateStr; /** \brief Pointer to the time generator. */ //GLG4TimeGen *fTimeGen; /** \brief Name of the time generator. */ // FIXME: Finish implementing ES support. For now this code is commented out. //std::map fVertexGenMap; // This is still causing some pain since to make it generic, the // ES generator also needs to be changed. For now keep just the IBD generator. // NeutrinoVertexGen *fVertexGen; /// Have a map of vertex generators (1 per spectrum type) /// This allows to work for a single spectrum type, and for multiple std::map fVertexGenMap; // Aux pointer to the generator that is currently being used to generate the vertex VertexGen_IBD * fCurrentVertexGen; /** \brief Position generator. Can be changed manually */ //GLG4PosGen *fPosGen; /** \brief Name of the low level generators. */ G4String fTimeGenName; G4String fPosGenName; G4String fVertexGenName; /** \brief Messenger for the solar generator */ ///FIXME: Check if it is actually needed // ReactorGenMessenger *fMessenger; private: // double NumNeutrinosPSecTotal(unsigned int ireactor); double NumNeutrinosPSecCore(unsigned int ireactor,unsigned int icore); static string PrintStateHelp(); string fVertexStateString; /// Mandatory parameters for the vertex generator /** \brief Flux to be generated. Possible options : "pp,pep,hep,be7,b8,c13,o15,f17". */ G4String fReactorListName; std::vector fReactorData; std::map fReactorSpectra; TH2D* fPowerPDF; // IBD is only sensitive to nu_e // ES is sensitive to both. Default to nue and ignore parameter if generator is IBD /** \brief Neutrino flavor to be generated. Possible options : "nue,numu,nuebar,numubar". * * This option is passed down to the vertex generator, which on the other hand passes it to the cross-section calculation. * \attention Keep in mind that 'numu' stands for the ES scattering cross-section for non-electron neutrinos (both muon and tau). * */ G4String fNuFlavor; /** \brief Flag to for using the ephemeris data. The properties have to be passed by messenger. */ G4bool fUseReactorPosition; /** \brief Direction of next incoming neutrino.*/ G4ThreeVector fNuDir; // Direction of the incoming neutrino for the next event to be generated /** EventTime object to deal with the non-sequentiality of the generated events. */ EventTime fEvTime; /** Flux scale in N SSMs*/ G4double fFluxScale; /** Event rate. It is the product between the flux and the total cross section * Needs to account for the flux scale. */ G4double fEvtRate; G4bool volumeInfoLoaded; /** Rate per target **/ G4double fRatePerTarget; /** Bool to apply reactor core power scaling as defined in REACTOR_STATUS tables */ G4bool applyCorePowerscaling; /// -- Auxiliary members G4String fSelectedReactor; G4int fSelectedCore; // This is hardcoded for now...but can safely be changed to a // higher value if there is a reactor with more than this number of cores static const int max_cores = 12; #ifdef REACTORDEV_DEBUG G4double fUTtime; #endif }; inline G4double ReactorGen::GetRandomNumber(const double rmin, const double rmax) { G4double rnd = G4UniformRand(); // random number from 0 to 1. G4double value = rmin + (rmax - rmin) * rnd; return value; } } #endif /* ReactorGen_H_ */