//////////////////////////////////////////////////////////////////// /// \class RAT::ESCrossSec /// /// \brief Calculates neutrino-electron elastic scattering. /// (based on original QSNO code by F. Duncan, M. Chen and Y. Takeuchi). /// /// \author Nuno Barros -- contact person /// /// REVISION HISTORY:\n /// 23-NOV-2010 - Nuno Barros - Imported and adapted into SNO+ RAT /// - Shamelessly taken from SNO QPhysics class PNuE /// /// 22-JUN-2012 - Nuno Barros - Cleaned up code and solved a small problem with unit conversion. /// /// 07-JUL-2012 - Nuno Barros - Removed all Geant4 streamers replacing them by RAT log objects. /// - Changed geant4 data types by C++ data types (suggestion by CIC). /// - Revised constness of members. /// /// 11-Oct-2019 - Nuno Barros /// * Modified code to maintain consistency with the CrossSectionBase class /// * Refactored code to use the cross section manager (maintains backwards compatibility) /// /// /// /// \details Calculates the neutrino-electron ES cross section. /// The calculation is performed either for the ES_e (W+Z) channel, or the ES_mu,tau (Z) channel. /// Some remarks concerning the calculations. /// \param E,Enu : Neutrino energy (MeV). /// \param T,Te : Recoil electron energy (MeV). /// \return Cross section in units of $10^{-42}cm^{2}$ /// /// Available strategies for the ES cross section calculation (set by messenger): /// - 1 : Original routine from QSNO::PNuE (Bahcall). /// - 2 : Improved routine from QSNO::PNuE (without rad. corrections). /// - 3 : Improved routine from QSNO::PNuE (with rad. corrections - analytical). /// - 4 (default) : Improved routine from QSNO::PNuE (with rad. corrections - table). /// /// //////////////////////////////////////////////////////////////////// #ifndef __RAT_ESCrossSec__ #define __RAT_ESCrossSec__ // G4 headers #include // RAT headers #include #include #include // ROOT headers #include namespace RAT { // Forward declarations within the namespace // class ESCrossSecMessenger; class ESCrossSec : public CrossSectionBase { public: enum NuEType {nue,nuebar,numu,numubar}; // Since we have multiple reactions, we cannot have static object declaration // Make these private to not use multiple instances // These should be private, right? // We would like to avoid having multiple ESCrossSec(const char* flavor = "nue"); virtual ~ESCrossSec(); /** * @brief Calculate the total cross section for the neutrino energy Enu. * * @param Enu * @return total cross section in units of \f$ 10^{-42} cm^{2} \f$ . */ virtual double Sigma(const double E) const; /** * @brief Calculate the differential cross section for the neutrino energy Enu. * * @param Enu Incoming neutrino energy (MeV). * @param Te Recoil electron energy (MeV). * @return Differencial cross section \f$ \frac{d\sigma}{dT} \f$ in units of \f$ 10^{-42} cm^{2} \f$ . */ double dSigmadT(const double Enu,const double Te) const; /** * Integrate the differential cross section between recoil energies T1 and T2. * * @param Enu Incoming neutrino energy (MeV). * @param T1 Lower limit of recoil energy interval (MeV). * @param T2 Upper limit of recoil energy interval (MeV). * @return \f$ \left.\frac{d\sigma}{dT}\right|_{\left[T1;T2\right]} \f$ in units of \f$ 10^{-42} cm^{2} \f$ . */ double IntegraldSigmadT(const double Enu,const double T1,const double T2) const ; /** * @brief Calculate the differential cross section as a function of the recoil angle. * * @param Enu Incoming neutrino energy. * @param CosTh Cosine of laboratory recoil angle of the electron. * @return \f$ \frac{d\sigma}{d\cos \theta} \f$ in units of \f$ 10^{-42} cm^{2} \f$ . */ double dSigmadCosTh(const double Enu,const double CosTh) const; /** * @brief 3D equivalent of ESCrossSec::dSigmadCosTh * * @param Enu Incoming neutrino energy. * @param theta Cosine (FIXME) of laboratory recoil angle of the electron. * @param phi Azimuthal recoil angle. * @return \f$ \frac{d\sigma}{d\Omega} \f$ in units of \f$ 10^{-42} cm^{2} \f$ . * * @note Usually not used. Needs further verification. * */ double dSigmadOmega(const double Enu,const double theta, const double phi) const; //----------------------------------------------------------------------- /** * @brief Setter for which calculation to use. * * The available calculations are: * -# 1 : No radiative corrections (Bahcall first calculation). * -# 2 : No radiative corrections from table. * -# 3 : Full analytical calculation with radiative corrections. * -# 4 : Full calculations with radiative corrections from table (default). * * @param strat Option for calculation to use. * */ virtual void SetStrategy(const short int strat); /** * @brief Getter of the cross section calculation type being used. * This is a wrapper to old functions, so that backwards compatility is * maintained. * * @return index of calculation type. */ inline short int GetRadiativeCorrection() const {return GetStrategy();} /** * @brief Setter for the Weak mixing angle. * * @param sintw Weak mixing angle : \f$ \sin \theta_{W} \f$ * */ void SetSinThetaW(const double &sintw); /** * @brief Getter for the Weak mixing angle. * * @return Weak mixing angle : \f$ \sin \theta_{W} \f$ * */ double GetSinThetaW() const {return fsinthetaW2;}; /** * @brief Setter for the neutrino type to be used. * * @param reaction neutrino type (nue,numu,nuebar,numubar) */ void SetReaction(const std::string &reaction); /** * @brief Getter for the neutrino type to be used. * * @return reaction neutrino type (nue,numu,nuebar,numubar) */ const std::string& GetReaction() const {return fReactionStr;}; // Override of base classes to set parameters in a RAT compatible way // This is used to call : // - SetSinThetaW // virtual void SetDPar(const std::string name, const double val); // This is used to call: // - SetStrategy (alternatively) virtual void SetIPar(const std::string name, const int val); /** * @brief Return a TGraph with the differential cross section for an incoming neutrino with energy Enu. * @param Enu Incoming neutrino energy (MeV) * @return TGraph with the shape of \f$ \frac{d\sigma}{dT} \f$ in units of \f$ 10^{-42} cm^{2} \f$ . */ TGraph *DrawdSigmadT(const double Enu) const; /** * Returns the global normalization of the cross section calculation. * For precision reasons, the cross-section is performed on a different scale, and therefore * any result returned by the calculation is missing this scale, which has to be applied separately. * * @return cross section scaling factor (1e-42) */ double CrossSecNorm() const {return 1e-42;}; protected: /** * @brief Internal function to load the data from the DB. */ virtual void LoadTablesDB(); void CalcG(); private: // Sets the defaults for the calculation void Defaults(); // //----------------------------------------------------------------------- // // add by y.t. 16-JAN-2003 double fL(const double x) const; // internal routine NuEType fReaction; /// Reaction type std::string fReactionStr; /// String characterizing the reaction type double fEmin,fEmax; /// Auxiliary variables to deal with the energies // Some constants static const double fGf; /// Fermi constant (GeV^-2) static const double fhbarc; /// hbar*c (MeV*fm) static const double fhbarc2; /// hbar*c^2(GeV^2 mb) static const double falpha; /// radiative correction term /** * This variable is defined as static (and not const) because it can be changed * in the macro file. However, the change should propagate to all instances of the class */ static double fsinthetaW2; double fMe; /// electron mass double fSigmaoverme; /// \f$ \frac{\sigma}{m_{e^{-}}}\f$ double fgL; double fgR; /// Vars to manipulate the tables. Since the tables are the same for all instances, these can be static static const double fTableTeMin; static const double fTableTeMax; // Total cross section table bool fExistTableTot; int fNDataTot; double fEnuStepTot; LinearInterp fTableTot_e; // Differential cross section table bool fExistTableDif; int fNDataDif; double fEnuStepDif; LinearInterp fTableDif_e; std::vector fTableDif; }; // // Some inline definitions to make the code a bit faster // inline void ESCrossSec::CalcG() { // calculate the gL and gR for // the reaction type if (fReaction == nue) { fgL = 0.5 + fsinthetaW2; fgR = fsinthetaW2; } else if (fReaction == nuebar) { fgL = fsinthetaW2; fgR = 0.5 + fsinthetaW2; } else if (fReaction == numu) { fgL = -0.5 + fsinthetaW2; fgR = fsinthetaW2; } else if (fReaction == numubar) { fgL = fsinthetaW2; fgR = -0.5 + fsinthetaW2; } else { throw std::invalid_argument("Unknown reaction " + fReactionStr); } } // //FIXME: This is not prepared to deal with multiple instances // // this is needed to simulate multiple flavors // inline ESCrossSec* // ESCrossSec::Get(const std::string& nu) // { // static ESCrossSec xs(nu.c_str()); // return &xs; // } // } // -- namespace RAT #endif