include/ninja/ninja.hh

0001 // Main -*- C++ -*- header file of the Ninja library.
0002
0003
0004 #ifndef NINJA_NINJA_HH
0005 #define NINJA_NINJA_HH
0006
0007 #define NINJA_NINJA_HH_INSIDE
0008
0009 #include <iostream>
0010 #include <stdexcept>
0011
0012 #include <ninja/types.hh>
0013 #include <ninja/momentum.hh>
0014 #include <ninja/s_mat.hh>
0015 #include <ninja/num_defs.hh>
0016 #include <ninja/ninja_in.hh>
0017
0018
0019 namespace ninja {
0020
0021   // Flags which control the internal tests performed by Ninja (NONE
0022   // by default).
0023   struct Test {
0024     enum {
0025       NONE = 0,
0026       ALL = ~(0),
0027       GLOBAL = 1,
0028       LOCAL_4 = 1 << 4,
0029       LOCAL_3 = 1 << 3,
0030       LOCAL_2 = 1 << 2,
0031       LOCAL_1 = 1 << 1,
0032       LOCAL = LOCAL_4 | LOCAL_3 | LOCAL_2 | LOCAL_1
0033     };
0034   };
0035
0036
0037   // Flags which control the output of ninja (NONE by default).
0038   struct Verbose {
0039     enum {
0040
0041       NONE = 0,
0042       ALL = ~(0),
0043
0044       // Output tests
0045       GLOBAL_TEST = 1,
0046       LOCAL_TEST_4 = 1 << 4,
0047       LOCAL_TEST_3 = 1 << 3,
0048       LOCAL_TEST_2 = 1 << 2,
0049       LOCAL_TEST_1 = 1 << 1,
0050       LOCAL_TESTS = LOCAL_TEST_1 | LOCAL_TEST_2 | LOCAL_TEST_3 | LOCAL_TEST_4,
0051       TESTS = LOCAL_TESTS | GLOBAL_TEST,
0052
0053       // Output coefficients
0054       C5 = 1 << (4+5),
0055       C4 = 1 << (4+4),
0056       C3 = 1 << (4+3),
0057       C2 = 1 << (4+2),
0058       C1 = 1 << (4+1),
0059       COEFFS = C5 | C4 | C3 | C2 | C1,
0060
0061       // Result
0062       RESULT = 1 << (4+5+1),
0063
0064       // Output integrals
0065       INTEGRALS = 1 << (4+5+2)
0066     };
0067   };
0068
0069   // Choose which tests to perform (NONE by default).  Example:
0070   // setTest(Test::LOCAL) performs local tests.
0071   inline void setTest(unsigned flag)
0072   {
0073     Options::test = flag;
0074   }
0075
0076   // Set the verbosity flag.  Example:
0077   // setVerbosity(Verbose::GLOBAL_TEST|Verbose::COEFFS) prints the
0078   // result of the global test and the value of all the computed
0079   // coefficients
0080   inline void setVerbosity(unsigned flag)
0081   {
0082     Options::verb = flag;
0083   }
0084
0085
0086   // set the integral library to be used by ninja
0087   inline void setDefaultIntegralLibrary(IntegralLibrary & mis)
0088   {
0089     Options::mis = & mis;
0090   }
0091
0092   // get a pointer to the integral library ninja is using
0093   inline IntegralLibrary * getIntegralLibrary()
0094   {
0095     return Options::mis;
0096   }
0097
0098   // Print the banner.  If force_print is set to false, the banner
0099   // will be printed only if never printed before.
0100   void printBanner(std::ostream & banner_out = std::cout,
0101                    bool force_print=true);
0102
0103   // Set the chop_tolerance.  If the verbosity is set to a non-zero
0104   // value, only the coefficients greater than chop_tolerance will be
0105   // printed as non-zeros.  Its default value is 1.0e-5.
0106   inline void setChopTolerance(Real chop_tolerance)
0107   {
0108     Options::chop_tol = chop_tolerance;
0109   }
0110
0111   // Set the chop_tolerance.  If a test is performed and the relative
0112   // error is greater than the specified tolerance, the
0113   // Amplitude::evaluate method will return Amplitude::TEST_FAILDED.
0114   // Its default value is 1.0e-5.
0115   inline void setTestTolerance(Real test_tolerance)
0116   {
0117     Options::test_tol = test_tolerance;
0118   }
0119
0120   // Set the default floating-point threshold used to detect unstable
0121   // kinematics to be used when not specified for an Amplitude object.
0122   inline void setDefaultFloatingPointThreshold(Real threshold)
0123   {
0124     Options::fp_threshold = threshold;
0125   }
0126
0127   // Set the output stream, used when the verbosity is true.
0128   inline void setOutputStream(std::ostream & outs)
0129   {
0130     Options::out = & outs;
0131   }
0132
0133
0134   // ninja::Amplitude is the main class of the ninja library.  The
0135   // method Amplitude::evaluate computes the integrals.
0136   template<typename MassType>
0137   class Amplitude {
0138   public:
0139
0140     enum {SUCCESS = 0, TEST_FAILED=1, UNSTABLE_KINEMATICS=1 << 1};
0141
0142     typedef typename const_pointer<MassType>::type MassConstPtr;
0143
0144   public:
0145
0146     // default
0147     Amplitude()
0148       : result(0.,0.,0.),
0149         s_mat(),
0150         V(0),
0151         mis(Options::mis),
0152         m2(static_cast<MassType*>(0)),
0153         cut_constr(),
0154         scale(1.),
0155         fp_threshold(Options::fp_threshold),
0156         n(0), rank(0),
0157         min_cut(0),
0158         return_val(SUCCESS),
0159         use_mu_exp(true) {}
0160
0161     // copy
0162     Amplitude(Amplitude & amp)
0163       : result(amp.result[0],amp.result[1],amp.result[2]),
0164         s_mat(),
0165         V(amp.V),
0166         mis(amp.mis),
0167         m2(amp.m2),
0168         cut_constr(amp.cut_constr),
0169         scale(amp.scale),
0170         fp_threshold(amp.fp_threshold),
0171         n(amp.n), rank(amp.rank),
0172         min_cut(amp.min_cut),
0173         return_val(SUCCESS),
0174         use_mu_exp(amp.use_mu_exp)
0175     {
0176       if (!amp.s_mat.isNull())
0177         s_mat.copy(amp.s_mat);
0178     }
0179
0180     // Constructor.  nn = number-of-loop-propagators,
0181     // rr=rank-of-numerator, p=pointer-to-array-ofinternal-momenta,
0182     // mass_sq=pointer-to-array-of-square-masses (which will be
0183     // ignored in the massless case))
0184     Amplitude(int nn, int rr, const RealMomentum p[], MassConstPtr mass_sq)
0185       : result(0.,0.,0.),
0186         s_mat(0,0),
0187         V(p),
0188         mis(Options::mis),
0189         m2(mass_sq),
0190         cut_constr(),
0191         scale(1.),
0192         fp_threshold(Options::fp_threshold),
0193         n(nn), rank(rr),
0194         min_cut(0),
0195         return_val(SUCCESS),
0196         use_mu_exp(true) {}
0197
0198     // this can only be used for MassType == Massless
0199     Amplitude(int nn, int rr, const RealMomentum p[]);
0200
0201     // Set number of loop propagators
0202     Amplitude & setN(int nn)
0203     {
0204       n = nn;
0205       return *this;
0206     }
0207
0208     // Set the rank of the numerator
0209     Amplitude & setRank(int rr)
0210     {
0211       rank = rr;
0212       return *this;
0213     }
0214
0215     // Set pointer to array of internal momenta
0216     Amplitude & setKinematics(const RealMomentum p[])
0217     {
0218       V = p;
0219       return *this;
0220     }
0221
0222     // Set pointer to array of internal squared masses
0223     Amplitude & setInternalMasses(MassConstPtr mass_sq)
0224     {
0225       m2 = mass_sq;
0226       return *this;
0227     }
0228
0229     // Set S-matrix from another
0230     Amplitude & setSMatrix(const SMatrix & s_matrix)
0231     {
0232       s_mat = s_matrix;
0233       return *this;
0234     }
0235
0236     // Set S-matrix from pointer to data
0237     Amplitude & setSMatrix(Real * data_ptr)
0238     {
0239       s_mat = SMatrix(data_ptr ? n : 0, data_ptr);
0240       return *this;
0241     }
0242
0243     // Set the renormalization scale (it only affects the computation
0244     // of the Master Integrals).  If not set, its default value is 1.
0245     Amplitude & setRenormalizationScale(Real renorm_scale)
0246     {
0247       scale = renorm_scale;
0248       return *this;
0249     }
0250
0251     // Stop the reduction right after the specified cut (e.g. if val=3
0252     // the reduction will stop after the determination of the triangle
0253     // residues and assumes that bubbles and tadpoles give no
0254     // contribution to the amplitude; if val<=1 instead the reduction
0255     // will be complete)
0256     Amplitude & setCutStop(int val)
0257     {
0258       min_cut = val;
0259       return *this;
0260     }
0261
0262     // Set the mu-expansion flag.  If set to true (default), use the
0263     // mu-expansion for the determination of the rational part coming
0264     // from the box contributions.  If set to false use the sampling
0265     // on mu^2 instead.
0266     Amplitude & useMuExpansion(bool val = true)
0267     {
0268       use_mu_exp = val;
0269       return *this;
0270     }
0271
0272     // set the integral library to be used by this Amplitude object,
0273     // if different from the default one
0274     Amplitude & setIntegralLibrary(IntegralLibrary & library)
0275     {
0276       mis = & library;
0277       return *this;
0278     }
0279
0280     // Reset the total result to zero
0281     Amplitude & reset()
0282     {
0283       result[0] = result[1] = result[2] = ninja::Complex();
0284       cut_constr = ninja::Complex();
0285       return *this;
0286     }
0287
0288     // Evaluate the Amplitude an add the result to the total.  Returns
0289     // Amplitude::SUCCESS unless one of the tests performed failed, in
0290     // which case it returns Amplitude::TEST_FAILED.
0291     int evaluate(Numerator & num);
0292
0293     // Acces to components eps^0, eps^-1, eps^-2 of the result
0294     const Complex eps0() const { return result[0];}
0295     const Complex epsm1() const { return result[1];}
0296     const Complex epsm2() const { return result[2];}
0297
0298     // rational part of the result
0299     const Complex getRationalPart() const
0300     {
0301       return result[0] - cut_constr;
0302     }
0303
0304     // cut constructible part of the result (finite-term)
0305     const Complex getCutConstructiblePart() const
0306     {
0307       return cut_constr;
0308     }
0309
0310     // throw away rational part
0311     void onlyCutConstructible()
0312     {
0313       result[0] = cut_constr;
0314     }
0315
0316     // operator []: operator [i] returns the coefficients of eps^(-i)
0317     const Complex operator[] (unsigned i) const
0318     {
0319       return result[i];
0320     }
0321
0322     // set floating-point threshold used to detect unstable
0323     // kinematics.  If an unstable kinematics is detected in the
0324     // 'evaluate' method, UNSTABLE_KINEMATICS is returned and the
0325     // result is left unchanged.
0326     void setFloatingPointThreshold(Real threshold)
0327     {
0328       fp_threshold = threshold;
0329     }
0330
0331   private:
0332
0333     void evaluatePentagons(Numerator & num,
0334                            CutsVector<cuts::Pentagon> & pentagon);
0335     void evaluatePentagon(Numerator & num,
0336                           cuts::Pentagon & pentagon);
0337
0338     void evaluateBoxes(Numerator & num,
0339                        const CutsVector<cuts::Pentagon> & pentagon,
0340                        CutsVector<cuts::Box> & box);
0341     void evaluateBox(Numerator & num,
0342                      const CutsVector<cuts::Pentagon> & pentagon,
0343                      cuts::Box & box);
0344
0345     void evaluateFullBoxes(Numerator & num,
0346                            const CutsVector<cuts::Pentagon> & pentagon,
0347                            CutsVector<cuts::Box> & box);
0348     void evaluateFullBox(Numerator & num,
0349                          const CutsVector<cuts::Pentagon> & pentagon,
0350                          cuts::Box & box);
0351
0352     void evaluateTriangles(Numerator & num,
0353                            CutsVector<cuts::Triangle> & triangle);
0354     void evaluateTriangle(Numerator & num,
0355                           cuts::Triangle & triangle);
0356
0357     void evaluateBubbles(Numerator & num,
0358                          const CutsVector<cuts::Triangle> & triangle,
0359                          CutsVector<cuts::Bubble> & bubble);
0360     void evaluateBubble(Numerator & num,
0361                         const CutsVector<cuts::Triangle> & triangle,
0362                         cuts::Bubble & bubble);
0363
0364     void evaluateTadpoles(Numerator & num,
0365                           const CutsVector<cuts::Triangle> & triangle,
0366                           const CutsVector<cuts::Bubble> & bubble,
0367                           CutsVector<cuts::Tadpole> & tadpole);
0368     void evaluateTadpole(Numerator & num,
0369                          const CutsVector<cuts::Triangle> & triangle,
0370                          const CutsVector<cuts::Bubble> & bubble,
0371                          cuts::Tadpole & tadpole);
0372
0373     void evaluateFullTadpoles(Numerator & num,
0374                               const CutsVector<cuts::Triangle> & triangle,
0375                               const CutsVector<cuts::Bubble> & bubble,
0376                               CutsVector<cuts::Tadpole> & tadpole);
0377     void evaluateFullTadpole(Numerator & num,
0378                              const CutsVector<cuts::Triangle> & triangle,
0379                              const CutsVector<cuts::Bubble> & bubble,
0380                              cuts::Tadpole & tadpole);
0381
0382     // Global N = N test
0383     void NeqNtest(Numerator & num,
0384                  const CutsVector<cuts::Pentagon> & pentagon,
0385                  const CutsVector<cuts::Box> & box,
0386                  const CutsVector<cuts::Triangle> & triangle,
0387                  const CutsVector<cuts::Bubble> & bubble,
0388                  const CutsVector<cuts::Tadpole> & tadpole,
0389                  const ComplexMomentum & q, const Complex & muq);
0390     // Local N = N tests
0391     void local4NeqNtests(Numerator & num,
0392                         const CutsVector<cuts::Pentagon> & pentagon,
0393                         const CutsVector<cuts::Box> & box);
0394     void local3NeqNtests(Numerator & num,
0395                         const CutsVector<cuts::Pentagon> & pentagon,
0396                         const CutsVector<cuts::Box> & box,
0397                         const CutsVector<cuts::Triangle> & triangle);
0398     void local2NeqNtests(Numerator & num,
0399                         const CutsVector<cuts::Pentagon> & pentagon,
0400                         const CutsVector<cuts::Box> & box,
0401                         const CutsVector<cuts::Triangle> & triangle,
0402                         const CutsVector<cuts::Bubble> & bubble);
0403     void local1NeqNtests(Numerator & num,
0404                         const CutsVector<cuts::Pentagon> & pentagon,
0405                         const CutsVector<cuts::Box> & box,
0406                         const CutsVector<cuts::Triangle> & triangle,
0407                         const CutsVector<cuts::Bubble> & bubble,
0408                         const CutsVector<cuts::Tadpole> & tadpole);
0409
0410
0411     // Higher rank methods
0412
0413     int higherRankEvaluate(Numerator & num);
0414
0415     void evaluatePentagons(Numerator & num,
0416                            CutsVector<x1cuts::Pentagon> & pentagon);
0417     void evaluatePentagon(Numerator & num,
0418                           x1cuts::Pentagon & pentagon);
0419
0420     void evaluateBoxes(Numerator & num,
0421                        const CutsVector<x1cuts::Pentagon> & pentagon,
0422                        CutsVector<x1cuts::Box> & box);
0423     void evaluateBox(Numerator & num,
0424                      const CutsVector<x1cuts::Pentagon> & pentagon,
0425                      x1cuts::Box & box);
0426
0427     void evaluateFullBoxes(Numerator & num,
0428                            const CutsVector<x1cuts::Pentagon> & pentagon,
0429                            CutsVector<x1cuts::Box> & box);
0430     void evaluateFullBox(Numerator & num,
0431                          const CutsVector<x1cuts::Pentagon> & pentagon,
0432                          x1cuts::Box & box);
0433
0434     void evaluateTriangles(Numerator & num,
0435                            CutsVector<x1cuts::Triangle> & triangle);
0436     void evaluateTriangle(Numerator & num,
0437                           x1cuts::Triangle & triangle);
0438
0439     void evaluateBubbles(Numerator & num,
0440                          const CutsVector<x1cuts::Triangle> & triangle,
0441                          CutsVector<x1cuts::Bubble> & bubble);
0442     void evaluateBubble(Numerator & num,
0443                         const CutsVector<x1cuts::Triangle> & triangle,
0444                         x1cuts::Bubble & bubble);
0445
0446     void evaluateTadpoles(Numerator & num,
0447                           const CutsVector<x1cuts::Triangle> & triangle,
0448                           const CutsVector<x1cuts::Bubble> & bubble,
0449                           CutsVector<x1cuts::Tadpole> & tadpole);
0450     void evaluateTadpole(Numerator & num,
0451                          const CutsVector<x1cuts::Triangle> & triangle,
0452                          const CutsVector<x1cuts::Bubble> & bubble,
0453                          x1cuts::Tadpole & tadpole);
0454
0455     void evaluateFullTadpoles(Numerator & num,
0456                               const CutsVector<x1cuts::Triangle> & triangle,
0457                               const CutsVector<x1cuts::Bubble> & bubble,
0458                               CutsVector<x1cuts::Tadpole> & tadpole);
0459     void evaluateFullTadpole(Numerator & num,
0460                              const CutsVector<x1cuts::Triangle> & triangle,
0461                              const CutsVector<x1cuts::Bubble> & bubble,
0462                              x1cuts::Tadpole & tadpole);
0463
0464     // If the argument is too small, set the kinematics as unstable
0465     // and returns false.  Returns true otherwise
0466     template <typename T>
0467     bool stability_check(const T & z)
0468     {
0469       if (taxicab_norm(z) < fp_threshold) {
0470         return_val = return_val | UNSTABLE_KINEMATICS;
0471         return false;
0472       }
0473       return true;
0474     }
0475
0476     // check if the kinematics has been detected as unstable
0477     bool unstable_kinematics() const
0478     {
0479       if (return_val & UNSTABLE_KINEMATICS) {
0480         if (Options::verb) {
0481           (*Options::out) << std::endl
0482                           << "ninja::Amplitude is returning "
0483                           << "UNSTABLE_KINEMATICS" << std::endl;
0484         }
0485         return true;
0486       }
0487       return false;
0488     }
0489
0490     // Global N = N test
0491     void NeqNtest(Numerator & num,
0492                  const CutsVector<x1cuts::Pentagon> & pentagon,
0493                  const CutsVector<x1cuts::Box> & box,
0494                  const CutsVector<x1cuts::Triangle> & triangle,
0495                  const CutsVector<x1cuts::Bubble> & bubble,
0496                  const CutsVector<x1cuts::Tadpole> & tadpole,
0497                  const ComplexMomentum & q, const Complex & muq);
0498     // Local N = N tests
0499     void local4NeqNtests(Numerator & num,
0500                         const CutsVector<x1cuts::Pentagon> & pentagon,
0501                         const CutsVector<x1cuts::Box> & box);
0502     void local3NeqNtests(Numerator & num,
0503                         const CutsVector<x1cuts::Pentagon> & pentagon,
0504                         const CutsVector<x1cuts::Box> & box,
0505                         const CutsVector<x1cuts::Triangle> & triangle);
0506     void local2NeqNtests(Numerator & num,
0507                         const CutsVector<x1cuts::Pentagon> & pentagon,
0508                         const CutsVector<x1cuts::Box> & box,
0509                         const CutsVector<x1cuts::Triangle> & triangle,
0510                         const CutsVector<x1cuts::Bubble> & bubble);
0511     void local1NeqNtests(Numerator & num,
0512                         const CutsVector<x1cuts::Pentagon> & pentagon,
0513                         const CutsVector<x1cuts::Box> & box,
0514                         const CutsVector<x1cuts::Triangle> & triangle,
0515                         const CutsVector<x1cuts::Bubble> & bubble,
0516                         const CutsVector<x1cuts::Tadpole> & tadpole);
0517
0518     // data
0519     ninja::details::Array3D<Complex> result;
0520     SMatrix s_mat;
0521     const RealMomentum * V;
0522     IntegralLibrary * mis;
0523     MassConstPtr m2;
0524     Complex cut_constr;
0525     Real scale, fp_threshold;
0526     int n, rank;
0527     int min_cut;
0528     int return_val;
0529     bool use_mu_exp;
0530
0531   }; // template<typename MassType> class Amplitude
0532
0533
0534   template <typename MassType>
0535   inline Amplitude<MassType>::Amplitude(int nn, int rr, const RealMomentum p[])
0536     : result(0.,0.,0.),
0537       s_mat(0,0),
0538       V(p),
0539       mis(Options::mis),
0540       m2(static_cast<const MassType*>(0)),
0541       cut_constr(),
0542       scale(1.),
0543       fp_threshold(Options::fp_threshold),
0544       n(nn), rank(rr),
0545       min_cut(0),
0546       return_val(SUCCESS),
0547       use_mu_exp(true)
0548   {
0549     // this makes sure it is only called for Massless types
0550     details::MasslessTypeError<MassType>::MassType_must_be_massless();
0551   }
0552
0553 } // namespace ninja
0554
0555 #undef NINJA_NINJA_HH_INSIDE
0556
0557 #endif // NINJA_NINJA_HH