//////////////////////////////////////////////////////////////////// /// \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. /// /// /// \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 // ROOT headers #include namespace RAT { // Forward declarations within the namespace class ESCrossSecMessenger; class ESCrossSec { public: enum NuEType {nue,nuebar,numu,numubar}; ESCrossSec(const char* flavor = "nue"); ~ESCrossSec(); // Set's the defaults for the calculation void Defaults(); /** * @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$ . */ double Sigma(const double Enu) 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 ii Option for calculation to use. * */ void SetRadiativeCorrection(const int ii); /** * @brief Getter of the cross section calculation type being used. * * @return index of calculation type. */ inline int GetRadiativeCorrection( ) const {return fRadiativeCorrection;}; /** * @brief Setter for the Weak mixing angle. * * @param sintw Weak mixing angle : \f$ \sin \theta_{W} \f$ * */ void SetSinThetaW(const double &sintw ) ; //{fsinthetaW2 = 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;}; /** * @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. */ void LoadTablesDB(); void CalcG(); private: /** * Rename of ESCrossSec::Sigma. To be discontinued. * * @see ESCrossSec::Sigma * @param Enu Incoming neutrino energy. * @return total cross section in units of \f$ 10^{-42} cm^{2} \f$ . */ double SigmaLab(double Enu) const; 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 sigmaoverme; /// \f$ \frac{\sigma}{m_{e^{-}}}\f$ double fgL; double fgR; //----------------------------------------------------------------------- // add by y.t. 16-JAN-2003 int fRadiativeCorrection; // flag to use radiative correction or not //double fL(double x); // internal routine /// 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 int fNDataTot; float fEnuStepTot; LinearInterp fTableTot_e; // Differential cross section table int fNDataDif; float fEnuStepDif; LinearInterp fTableDif_e; std::vector fTableDif; ESCrossSecMessenger *fMessenger; }; // // 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); } } #endif