AMISIC%2B%2B/Tools/Interaction_Probability.H

0001 #ifndef AMISIC_Tools_Interaction_Probability_H
0002 #define AMISIC_Tools_Interaction_Probability_H
0003
0004 #include "AMISIC++/Tools/Matter_Overlap.H"
0005 #include "AMISIC++/Tools/Lookup_Tables.H"
0006
0007 namespace AMISIC {
0008   class MI_Processes;
0009
0010   class Interaction_Probability {
0011     Matter_Overlap   m_mo;
0012     OneDim_Table   * p_k, * p_integral, * p_expO, * p_fc;
0013     bool   m_test;
0014
0015     void   FixK(MI_Processes * processes);
0016     double Integral(const double & k,const int & diff=0);
0017     double NewtonRaphson(const double & xsratio);
0018     void   FixOExp();
0019     void   OutputTables(MI_Processes * processes);
0020   public:
0021     Interaction_Probability();
0022     ~Interaction_Probability();
0023
0024     void   Initialize(REMNANTS::Remnant_Handler * remnant_handler,
0025               MI_Processes * processes);
0026
0027     inline Matter_Overlap * GetOverlap()                        { return &m_mo; }
0028     /////////////////////////////////////////////////////////////////////////////////
0029     // Probability for at least one interaction, Eq. (24)
0030     /////////////////////////////////////////////////////////////////////////////////
0031     inline double operator()(const double & s,const double & b) { return 1.-exp(-(*p_k)(s) * m_mo(b)); }
0032     /////////////////////////////////////////////////////////////////////////////////
0033     // b-dependent modification factor, Eq. (28)
0034     /////////////////////////////////////////////////////////////////////////////////
0035     inline double fb(const double & s,const double & b)         { return m_mo(b)/(*p_expO)(s); }
0036     /////////////////////////////////////////////////////////////////////////////////
0037     // b-dependent modification factor, Eq. (28)
0038     /////////////////////////////////////////////////////////////////////////////////
0039     inline double expO(const double & s)                        { return (*p_expO)(s); }
0040     /////////////////////////////////////////////////////////////////////////////////
0041     // Correction factor for demaning at least one interaction, Eq. (31)
0042     /////////////////////////////////////////////////////////////////////////////////
0043     inline double fc(const double & s)                          { return (*p_fc)(s); }
0044   };
0045
0046
0047   /////////////////////////////////////////////////////////////////////////////////
0048   // Interaction probability, Eq. (24) as operator, needs to be integrated for
0049   // Eqs.(26, 32)
0050   /////////////////////////////////////////////////////////////////////////////////
0051   class P_Integrand : public ATOOLS::Function_Base {
0052     Matter_Overlap * p_mo;
0053     double           m_k;
0054   public:
0055     P_Integrand(Matter_Overlap * mo, const double & k) : p_mo(mo), m_k(k) {}
0056     ~P_Integrand() {};
0057     double operator()(double b);
0058   };
0059
0060   /////////////////////////////////////////////////////////////////////////////////
0061   // Integrand O(b) exp[-kO(b)], used by the Newton-Raphson method to fix k
0062   /////////////////////////////////////////////////////////////////////////////////
0063   class OtimesExp_Integrand : public ATOOLS::Function_Base {
0064     Matter_Overlap * p_mo;
0065     double           m_k;
0066   public:
0067     OtimesExp_Integrand(Matter_Overlap * mo, const double & k) : p_mo(mo), m_k(k) {}
0068     ~OtimesExp_Integrand() {};
0069     double operator()(double b);
0070   };
0071
0072   /////////////////////////////////////////////////////////////////////////////////
0073   // Integrand of numerator in Eqs. (29, 31) as operator, given by O(b) Pint(b)
0074   /////////////////////////////////////////////////////////////////////////////////
0075   class OtimesP_Integrand : public ATOOLS::Function_Base {
0076     Matter_Overlap * p_mo;
0077     double           m_k;
0078   public:
0079     OtimesP_Integrand(Matter_Overlap * mo, const double & k) : p_mo(mo), m_k(k) {}
0080     ~OtimesP_Integrand() {};
0081     double operator()(double b);
0082   };
0083 }
0084
0085 #endif