Matchbox/Phasespace/RandomHelpers.h

0001 // -*- C++ -*-
0002 //
0003 // RandomHelpers.h is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef HERWIG_RandomHelpers_H
0010 #define HERWIG_RandomHelpers_H
0011
0012 #include "ThePEG/Config/ThePEG.h"
0013
0014 namespace Herwig {
0015
0016 using namespace ThePEG;
0017
0018 /**
0019  * \ingroup Matchbox
0020  * \author Simon Platzer
0021  * \brief Phase space generation utilities.
0022  */
0023 namespace RandomHelpers {
0024
0025 /**
0026  * \ingroup Matchbox
0027  * \author Simon Platzer
0028  * \brief Small helper.
0029  */
0030 inline double sign(double x) {
0031   return x < 0. ? -1. : 1.;
0032 }
0033
0034 /**
0035  * \ingroup Matchbox
0036  * \author Simon Platzer
0037  * \brief Define the generator concept.
0038  */
0039 template<class Density>
0040 struct Generator {
0041
0042   /**
0043    * Return the lower bound of the density generated.
0044    */
0045   double lower() const;
0046
0047   /**
0048    * Return the upper bound of the density generated.
0049    */
0050   double upper() const;
0051
0052   /**
0053    * Return the density's value
0054    */
0055   double value(double x) const;
0056
0057   /**
0058    * Return the density's normalization
0059    */
0060   double normalization() const;
0061
0062   /**
0063    * Generate the return value according to the implemented density,
0064    * given a flat random number on the unit interval.
0065    */
0066   double operator()(double r) const;
0067
0068 };
0069
0070 /**
0071  * \ingroup Matchbox
0072  * \author Simon Platzer
0073  * \brief A density expression.
0074  */
0075 struct Expression {};
0076
0077 /**
0078  * \ingroup Matchbox
0079  * \author Simon Platzer
0080  * \brief Container base class for a general density.
0081  */
0082 template<>
0083 struct Generator<Expression> {
0084
0085   /**
0086    * The destructor.
0087    */
0088   virtual ~Generator() {}
0089
0090   /**
0091    * Return the lower bound of the density generated.
0092    */
0093   virtual double lower() const = 0;
0094
0095   /**
0096    * Return the upper bound of the density generated.
0097    */
0098   virtual double upper() const = 0;
0099
0100   /**
0101    * Return the density's value
0102    */
0103   virtual double value(double x) const = 0;
0104
0105   /**
0106    * Return the density's normalization
0107    */
0108   virtual double normalization() const = 0;
0109
0110   /**
0111    * Generate the return value according to the implemented density,
0112    * given a flat random number on the unit interval.
0113    */
0114   virtual double operator()(double r) const = 0;
0115
0116 };
0117
0118 /**
0119  * \ingroup Matchbox
0120  * \author Simon Platzer
0121  * \brief A density container.
0122  */
0123 template<class Density>
0124 struct Container {};
0125
0126 /**
0127  * \ingroup Matchbox
0128  * \author Simon Platzer
0129  * \brief Container class for a general density.
0130  */
0131 template<class Density>
0132 class Generator<Container<Density> >
0133   : public Generator<Expression> {
0134
0135   /**
0136    * The generator.
0137    */
0138   Generator<Density> generator;
0139
0140 public:
0141
0142   /**
0143    * Construct from generator.
0144    */
0145   Generator(const Generator<Density>& gen)
0146     : generator(gen) {}
0147
0148   /**
0149    * Return the lower bound of the density generated.
0150    */
0151   virtual double lower() const { return generator.lower(); }
0152
0153   /**
0154    * Return the upper bound of the density generated.
0155    */
0156   virtual double upper() const { return generator.upper(); }
0157
0158   /**
0159    * Return the density's value
0160    */
0161   virtual double value(double x) const { return generator.value(x); }
0162
0163   /**
0164    * Return the density's normalization
0165    */
0166   virtual double normalization() const { return generator.normalization(); }
0167
0168   /**
0169    * Generate the return value according to the implemented density,
0170    * given a flat random number on the unit interval.
0171    */
0172   virtual double operator()(double r) const { return generator(r); }
0173
0174 };
0175
0176 /**
0177  * \ingroup Matchbox
0178  * \author Simon Platzer
0179  * \brief Generate a random variable and return its weight.
0180  */
0181 template<class Density>
0182 pair<double,double> generate(const Generator<Density>& gen,
0183                  double r) {
0184   double x = gen(r);
0185
0186   if ( gen.value(x) != 0. )
0187     return make_pair(x,gen.normalization()/gen.value(x));
0188   else
0189     return make_pair(x,0.);
0190 }
0191
0192 /**
0193  * \ingroup Matchbox
0194  * \author Simon Platzer
0195  * \brief Remap a density to a new interval.
0196  */
0197 template<class Density>
0198 struct Remap {};
0199
0200 /**
0201  * \ingroup Matchbox
0202  * \author Simon Platzer
0203  * \brief Generate a density remapped to a new interval.
0204  */
0205 template<class Density>
0206 class Generator<Remap<Density> > {
0207
0208   /**
0209    * The underlying generator.
0210    */
0211   Generator<Density> theGenerator;
0212
0213   /**
0214    * The new lower bound.
0215    */
0216   double theLower;
0217
0218   /**
0219    * The new upper bound.
0220    */
0221   double theUpper;
0222
0223   /**
0224    * Construct from generator and new boundaries.
0225    */
0226   Generator(const Generator<Density>& gen,
0227         double low, double up)
0228     : theGenerator(gen), theLower(low), theUpper(up) {
0229     if ( low >= up )
0230       throw std::logic_error("[Generator<Remap>] Invalid boundaries.");
0231   }
0232
0233   /**
0234    * Return the generator.
0235    */
0236   const Generator<Density>& generator() const { return theGenerator; }
0237
0238   /**
0239    * Return the lower bound of the density generated.
0240    */
0241   double lower() const { return theLower; }
0242
0243   /**
0244    * Return the upper bound of the density generated.
0245    */
0246   double upper() const { return theUpper; }
0247
0248   /**
0249    * Return the density's value
0250    */
0251   double value(double y) const {
0252     double xm = generator().lower();
0253     double xp = generator().upper();
0254     double ym = lower();
0255     double yp = upper();
0256     double x = ((xp-xm)/(yp-ym))*y+(yp*xm-ym*xp)/(yp-ym);
0257     return generator().value(x);
0258   }
0259
0260   /**
0261    * Return the density's normalization
0262    */
0263   double normalization() const {
0264     double xm = generator().lower();
0265     double xp = generator().upper();
0266     double ym = lower();
0267     double yp = upper();
0268     return ((yp-ym)/(xp-xm))*generator().normalization();
0269   }
0270
0271   /**
0272    * Generate the return value according to the implemented density,
0273    * given a flat random number on the unit interval.
0274    */
0275   double operator()(double r) const {
0276     double xm = generator().lower();
0277     double xp = generator().upper();
0278     double ym = lower();
0279     double yp = upper();
0280     double x = ((yp-ym)/(xp-xm))*generator()(r)+(xp*ym-xm*yp)/(xp-xm);
0281     return x;
0282   }
0283
0284 };
0285
0286 /**
0287  * \ingroup Matchbox
0288  * \author Simon Platzer
0289  * \brief Indicate remapping of a density.
0290  */
0291 struct on {
0292
0293   /**
0294    * The new lower boundary.
0295    */
0296   double lower;
0297
0298   /**
0299    * The new upper boundary.
0300    */
0301   double upper;
0302
0303   /**
0304    * Construct from boundaries.
0305    */
0306   on(double a, double b)
0307     :lower(a), upper(b) {}
0308 };
0309
0310 /**
0311  * \ingroup Matchbox
0312  * \author Simon Platzer
0313  * \brief Construct a remapped density generator.
0314  */
0315 template<class Density>
0316 Generator<Remap<Density> > operator*(const Generator<Density>& gen,
0317                      const on& interval) {
0318   return Generator<Remap<Density> >(gen,interval.lower,interval.upper);
0319 }
0320
0321
0322 /**
0323  * \ingroup Matchbox
0324  * \author Simon Platzer
0325  * \brief Rescale a density.
0326  */
0327 template<class Density>
0328 struct Rescale {};
0329
0330 /**
0331  * \ingroup Matchbox
0332  * \author Simon Platzer
0333  * \brief Generate a rescaled density.
0334  */
0335 template<class Density>
0336 class Generator<Rescale<Density> > {
0337
0338   /**
0339    * The underlying generator.
0340    */
0341   Generator<Density> theGenerator;
0342
0343   /**
0344    * The rescaling factor.
0345    */
0346   double theScale;
0347
0348 public:
0349
0350   /**
0351    * Construct from generator and scale.
0352    */
0353   Generator(const Generator<Density>& gen,
0354         double sc)
0355     : theGenerator(gen), theScale(sc) {
0356   }
0357
0358   /**
0359    * Return the generator.
0360    */
0361   const Generator<Density>& generator() const { return theGenerator; }
0362
0363   /**
0364    * Return the scale
0365    */
0366   double scale() const { return theScale; }
0367
0368   /**
0369    * Return the lower bound of the density generated.
0370    */
0371   double lower() const { return generator().lower(); }
0372
0373   /**
0374    * Return the upper bound of the density generated.
0375    */
0376   double upper() const { return generator().upper(); }
0377
0378   /**
0379    * Return the density's value
0380    */
0381   double value(double x) const {
0382     return scale()*generator().value(x);
0383   }
0384
0385   /**
0386    * Return the density's normalization
0387    */
0388   double normalization() const {
0389     return scale()*generator().normalization();
0390   }
0391
0392   /**
0393    * Generate the return value according to the implemented density,
0394    * given a flat random number on the unit interval.
0395    */
0396   double operator()(double r) const {
0397     return generator()(r);
0398   }
0399
0400 };
0401
0402 /**
0403  * \ingroup Matchbox
0404  * \author Simon Platzer
0405  * \brief Construct a rescaled density.
0406  */
0407 template<class Density>
0408 Generator<Rescale<Density> > operator*(double a, const Generator<Density>& gen) {
0409   return Generator<Rescale<Density> >(gen,a);
0410 }
0411
0412 /**
0413  * \ingroup Matchbox
0414  * \author Simon Platzer
0415  * \brief Add two densities.
0416  */
0417 template<class Density1,
0418      class Density2>
0419 struct Sum {};
0420
0421 /**
0422  * \ingroup Matchbox
0423  * \author Simon Platzer
0424  * \brief Generate the sum of two densities.
0425  */
0426 template<class Density1,
0427      class Density2>
0428 class Generator<Sum<Density1,Density2> > {
0429
0430   /**
0431    * The first generator.
0432    */
0433   Generator<Density1> theFirstGenerator;
0434
0435   /**
0436    * The second generator.
0437    */
0438   Generator<Density2> theSecondGenerator;
0439
0440   /**
0441    * The lower boundary.
0442    */
0443   double theLower;
0444
0445   /**
0446    * The upper boundary.
0447    */
0448   double theUpper;
0449
0450   /**
0451    * The fraction of the unit interval considered for the first
0452    * generator.
0453    */
0454   double theFraction;
0455
0456 public:
0457
0458   /**
0459    * Construct from generators.
0460    */
0461   Generator(const Generator<Density1>& firstGen,
0462         const Generator<Density2>& secondGen)
0463     : theFirstGenerator(firstGen), theSecondGenerator(secondGen),
0464       theLower(min(firstGen.lower(),secondGen.lower())),
0465       theUpper(max(firstGen.upper(),secondGen.upper())),
0466       theFraction(1.) {
0467     theFraction =
0468       firstGenerator().normalization() / normalization();
0469   }
0470
0471   /**
0472    * Return the first generator.
0473    */
0474   const Generator<Density1>& firstGenerator() const { return theFirstGenerator; }
0475
0476   /**
0477    * Return the second generator.
0478    */
0479   const Generator<Density2>& secondGenerator() const { return theSecondGenerator; }
0480
0481   /**
0482    * Return the lower bound of the density generated.
0483    */
0484   double lower() const { return theLower; }
0485
0486   /**
0487    * Return the upper bound of the density generated.
0488    */
0489   double upper() const { return theUpper; }
0490
0491   /**
0492    * Return the fraction of the unit interval considered for the first
0493    * generator.
0494    */
0495   double fraction() const { return theFraction; }
0496
0497   /**
0498    * Return the density's value
0499    */
0500   double value(double x) const {
0501     double res = 0.;
0502     if ( firstGenerator().lower() <= x &&
0503      x <= firstGenerator().upper() )
0504       res += firstGenerator().value(x);
0505     if ( secondGenerator().lower() <= x &&
0506      x <= secondGenerator().upper() )
0507       res += secondGenerator().value(x);
0508     return res;
0509   }
0510
0511   /**
0512    * Return the density's normalization
0513    */
0514   double normalization() const {
0515     return
0516       firstGenerator().normalization() + secondGenerator().normalization();
0517   }
0518
0519   /**
0520    * Generate the return value according to the implemented density,
0521    * given a flat random number on the unit interval.
0522    */
0523   double operator()(double r) const {
0524     return
0525       r < fraction() ?
0526       firstGenerator()(r/fraction()) :
0527       secondGenerator()((r-fraction())/(1.-fraction()));
0528   }
0529
0530 };
0531
0532 /**
0533  * \ingroup Matchbox
0534  * \author Simon Platzer
0535  * \brief Construct the sum of two densities.
0536  */
0537 template<class Density1,
0538      class Density2>
0539 Generator<Sum<Density1,Density2> > operator+(const Generator<Density1>& first,
0540                          const Generator<Density2>& second) {
0541   return Generator<Sum<Density1,Density2> >(first,second);
0542 }
0543
0544 /**
0545  * \ingroup Matchbox
0546  * \author Simon Platzer
0547  * \brief Indicate that the argument density should be matched to the
0548  * previous one in a piecewise definition.
0549  */
0550 template<class Density>
0551 struct matcher {
0552
0553   /**
0554    * The generator to be matched.
0555    */
0556   Generator<Density> generator;
0557
0558   /**
0559    * Construct from generator.
0560    */
0561   matcher(const Generator<Density>& gen)
0562     : generator(gen) {}
0563
0564 };
0565
0566 /**
0567  * \ingroup Matchbox
0568  * \author Simon Platzer
0569  * \brief Indicate that the argument density should be matched to the
0570  * previous one in a piecewise definition.
0571  */
0572 template<class Density>
0573 matcher<Density> match(const Generator<Density>& gen) {
0574   return matcher<Density>(gen);
0575 }
0576
0577 /**
0578  * \ingroup Matchbox
0579  * \author Simon Platzer
0580  * \brief Construct the sum of two densities, matching the first
0581  * summand at its upper bound to the second at its lower bound.
0582  */
0583 template<class Density1,
0584      class Density2>
0585 Generator<Sum<Density1,Rescale<Density2> > > operator+(const Generator<Density1>& first,
0586                                const matcher<Density2>& second) {
0587   double matching =
0588     first.value(first.upper())/
0589     second.generator.value(second.generator.lower());
0590   return Generator<Sum<Density1,Rescale<Density2> > >(first,matching*second.generator);
0591 }
0592
0593 /**
0594  * \ingroup Matchbox
0595  * \author Simon Platzer
0596  * \brief A piecewise defined density.
0597  */
0598 template<class Density1,
0599      class Density2>
0600 struct Piecewise {};
0601
0602 /**
0603  * \ingroup Matchbox
0604  * \author Simon Platzer
0605  * \brief Placeholder when constructing piecewise defined densities.
0606  */
0607 struct ToBeDefined {};
0608
0609 /**
0610  * \ingroup Matchbox
0611  * \author Simon Platzer
0612  * \brief Generate a piecewise defined density.
0613  */
0614 template<class Density1,
0615      class Density2>
0616 class Generator<Piecewise<Density1,Density2> > {
0617
0618   /**
0619    * The first generator.
0620    */
0621   Generator<Density1> theFirstGenerator;
0622
0623   /**
0624    * The second generator.
0625    */
0626   Generator<Density2> theSecondGenerator;
0627
0628   /**
0629    * The lower boundary.
0630    */
0631   double theLower;
0632
0633   /**
0634    * The transition value.
0635    */
0636   double theIntermediate;
0637
0638   /**
0639    * The upper boundary.
0640    */
0641   double theUpper;
0642
0643   /**
0644    * The fraction of the unit interval considered for the first
0645    * generator.
0646    */
0647   double theFraction;
0648
0649 public:
0650
0651   /**
0652    * Construct from generators.
0653    */
0654   Generator(const Generator<Density1>& firstGen,
0655         const Generator<Density2>& secondGen)
0656     : theFirstGenerator(firstGen), theSecondGenerator(secondGen),
0657       theLower(firstGen.lower()), theIntermediate(firstGen.upper()), theUpper(secondGen.upper()),
0658       theFraction(1.) {
0659     if ( firstGenerator().upper() != secondGenerator().lower() )
0660       throw std::logic_error("[Generator<Piecewise>] Invalid boundaries.");
0661     theFraction =
0662       firstGenerator().normalization() / normalization();
0663   }
0664
0665   /**
0666    * Return the first generator.
0667    */
0668   const Generator<Density1>& firstGenerator() const { return theFirstGenerator; }
0669
0670   /**
0671    * Return the second generator.
0672    */
0673   const Generator<Density2>& secondGenerator() const { return theSecondGenerator; }
0674
0675   /**
0676    * Return the lower bound of the density generated.
0677    */
0678   double lower() const { return theLower; }
0679
0680   /**
0681    * Return the transition value.
0682    */
0683   double intermediate() const { return theIntermediate; }
0684
0685   /**
0686    * Return the upper bound of the density generated.
0687    */
0688   double upper() const { return theUpper; }
0689
0690   /**
0691    * Return the fraction of the unit interval considered for the first
0692    * generator.
0693    */
0694   double fraction() const { return theFraction; }
0695
0696   /**
0697    * Return the density's value
0698    */
0699   double value(double x) const {
0700     return
0701       x < intermediate() ?
0702       firstGenerator().value(x) :
0703       secondGenerator().value(x);
0704   }
0705
0706   /**
0707    * Return the density's normalization
0708    */
0709   double normalization() const {
0710     return
0711       firstGenerator().normalization() + secondGenerator().normalization();
0712   }
0713
0714   /**
0715    * Generate the return value according to the implemented density,
0716    * given a flat random number on the unit interval.
0717    */
0718   double operator()(double r) const {
0719     return
0720       r < fraction() ?
0721       firstGenerator()(r/fraction()) :
0722       secondGenerator()((r-fraction())/(1.-fraction()));
0723   }
0724
0725   /**
0726    * Construct piecewise generators.
0727    */
0728   template<class Density3>
0729   Generator<Piecewise<Piecewise<Density1,Density2>,Density3> >
0730   operator,(const Generator<Density3>& thirdGenerator) {
0731     return
0732       Generator<Piecewise<Piecewise<Density1,Density2>,Density3> >
0733       (*this,thirdGenerator);
0734   }
0735
0736   /**
0737    * Construct piecewise generators.
0738    */
0739   template<class Density3>
0740   Generator<Piecewise<Piecewise<Density1,Density2>,Rescale<Density3> > >
0741   operator,(const matcher<Density3>& thirdGenerator) {
0742     double matching =
0743       value(upper())/thirdGenerator.generator.value(upper());
0744     return
0745       Generator<Piecewise<Piecewise<Density1,Density2>,Rescale<Density3> > >
0746       (*this,matching*thirdGenerator.generator);
0747   }
0748
0749 };
0750
0751 /**
0752  * \ingroup Matchbox
0753  * \author Simon Platzer
0754  * \brief Generate a piecewise defined density.
0755  */
0756 template<class Density>
0757 struct Generator<Piecewise<Density,ToBeDefined> > {
0758
0759   /**
0760    * The first generator.
0761    */
0762   Generator<Density> generator;
0763
0764   /**
0765    * Construct from generator.
0766    */
0767   Generator(const Generator<Density>& gen)
0768     : generator(gen) {}
0769
0770   /**
0771    * Construct piecewise generators.
0772    */
0773   template<class Density2>
0774   Generator<Piecewise<Density,Density2> >
0775   operator,(const Generator<Density2>& secondGen) {
0776     return
0777       Generator<Piecewise<Density,Density2> >
0778       (generator,secondGen);
0779   }
0780
0781   /**
0782    * Construct piecewise generators.
0783    */
0784   template<class Density2>
0785   Generator<Piecewise<Density,Rescale<Density2> > >
0786   operator,(const matcher<Density2>& secondGen) {
0787     double matching =
0788       generator.value(generator.upper())/secondGen.generator.value(generator.upper());
0789     return
0790       Generator<Piecewise<Density,Rescale<Density2> > >
0791       (generator,matching*secondGen.generator);
0792   }
0793
0794 };
0795
0796 /**
0797  * \ingroup Matchbox
0798  * \author Simon Platzer
0799  * \brief Generate a piecewise defined density.
0800  */
0801 template<>
0802 struct Generator<Piecewise<ToBeDefined,ToBeDefined> > {
0803
0804   /**
0805    * Construct piecewise generators.
0806    */
0807   template<class Density>
0808   Generator<Piecewise<Density,ToBeDefined> >
0809   operator,(const Generator<Density>& gen) {
0810     return
0811       Generator<Piecewise<Density,ToBeDefined> >(gen);
0812   }
0813
0814 };
0815
0816 /**
0817  * \ingroup Matchbox
0818  * \author Simon Platzer
0819  * \brief Construct a piecewise defined density.
0820  */
0821 inline Generator<Piecewise<ToBeDefined,ToBeDefined> >
0822 piecewise() {
0823   return Generator<Piecewise<ToBeDefined,ToBeDefined> >();
0824 }
0825
0826
0827
0828 /**
0829  * \ingroup Matchbox
0830  * \author Simon Platzer
0831  * \brief A constant density.
0832  */
0833 struct Flat {};
0834
0835 /**
0836  * \ingroup Matchbox
0837  * \author Simon Platzer
0838  * \brief Generate x flat.
0839  */
0840 template<>
0841 class Generator<Flat> {
0842
0843   /**
0844    * The lower boundary.
0845    */
0846   double theLower;
0847
0848   /**
0849    * The upper boundary.
0850    */
0851   double theUpper;
0852
0853 public:
0854
0855   /**
0856    * Construct from boundaries.
0857    */
0858   Generator(double low, double up)
0859     : theLower(low), theUpper(up) {}
0860
0861   /**
0862    * Return the lower bound of the density generated.
0863    */
0864   double lower() const { return theLower; }
0865
0866   /**
0867    * Return the upper bound of the density generated.
0868    */
0869   double upper() const { return theUpper; }
0870
0871   /**
0872    * Return the density's value
0873    */
0874   double value(double x) const {
0875     return x>=lower() && x<=upper() ? 1. : 0.;
0876   }
0877
0878   /**
0879    * Return the density's normalization
0880    */
0881   double normalization() const { return upper()-lower(); }
0882
0883   /**
0884    * Generate the return value according to the implemented density,
0885    * given a flat random number on the unit interval.
0886    */
0887   double operator()(double r) const {
0888     return lower() + r*(upper()-lower());
0889   }
0890
0891 };
0892
0893 /**
0894  * \ingroup Matchbox
0895  * \author Simon Platzer
0896  * \brief Construct a constant density.
0897  */
0898 inline Generator<Flat> flat(double low, double up) {
0899   return Generator<Flat>(low,up);
0900 }
0901
0902 /**
0903  * \ingroup Matchbox
0904  * \author Simon Platzer
0905  * \brief A zero density.
0906  */
0907 struct Zero {};
0908
0909 /**
0910  * \ingroup Matchbox
0911  * \author Simon Platzer
0912  * \brief Generate nothing.
0913  */
0914 template<>
0915 class Generator<Zero> {
0916
0917   /**
0918    * The lower boundary.
0919    */
0920   double theLower;
0921
0922   /**
0923    * The upper boundary.
0924    */
0925   double theUpper;
0926
0927 public:
0928
0929   /**
0930    * Construct from boundaries.
0931    */
0932   Generator(double low, double up)
0933     : theLower(low), theUpper(up) {}
0934
0935   /**
0936    * Return the lower bound of the density generated.
0937    */
0938   double lower() const { return theLower; }
0939
0940   /**
0941    * Return the upper bound of the density generated.
0942    */
0943   double upper() const { return theUpper; }
0944
0945   /**
0946    * Return the density's value
0947    */
0948   double value(double x) const {
0949     return x>=lower() && x<=upper() ? Constants::epsilon : 0.;
0950   }
0951
0952   /**
0953    * Return the density's normalization
0954    */
0955   double normalization() const { return 0.; }
0956
0957   /**
0958    * Generate the return value according to the implemented density,
0959    * given a flat random number on the unit interval.
0960    */
0961   double operator()(double r) const {
0962     return lower() + r*(upper()-lower());
0963   }
0964
0965 };
0966
0967 /**
0968  * \ingroup Matchbox
0969  * \author Simon Platzer
0970  * \brief Construct a zero density.
0971  */
0972 inline Generator<Zero> zero(double low, double up) {
0973   return Generator<Zero>(low,up);
0974 }
0975
0976 /**
0977  * \ingroup Matchbox
0978  * \author Simon Platzer
0979  * \brief A density 1/|x-z|
0980  */
0981 struct Inverse {};
0982
0983 /**
0984  * \ingroup Matchbox
0985  * \author Simon Platzer
0986  * \brief Generate x with density 1/|x-z|
0987  */
0988 template<>
0989 class Generator<Inverse> {
0990
0991   /**
0992    * The position of the pole
0993    */
0994   double thePole;
0995
0996   /**
0997    * The lower bound
0998    */
0999   double theLower;
1000
1001   /**
1002    * The upper bound
1003    */
1004   double theUpper;
1005
1006   /**
1007    * Scale for random numbers.
1008    */
1009   double theScale;
1010
1011   /**
1012    * Offset for random mnumbers.
1013    */
1014   double theOffset;
1015
1016 public:
1017
1018   /**
1019    * Construct from pole and boundaries.
1020    */
1021   Generator(double z,
1022         double l, double u)
1023     : thePole(z),
1024       theLower(l), theUpper(u),
1025       theScale(z < l ? log((u-z)/(l-z)) : log((z-l)/(z-u))),
1026       theOffset(z < l ? log(l-z) : log(z-u)) {
1027     if ( z >= l && z <= u )
1028       throw std::logic_error("[Generator<Inverse>] Pole inside sampling interval.");
1029   }
1030
1031   /**
1032    * Return the lower bound of the density generated.
1033    */
1034   double lower() const { return theLower; }
1035
1036   /**
1037    * Return the upper bound of the density generated.
1038    */
1039   double upper() const { return theUpper; }
1040
1041   /**
1042    * Return the position of the pole.
1043    */
1044   double pole() const { return thePole; }
1045
1046   /**
1047    * Return the scale for random numbers.
1048    */
1049   double scale() const { return theScale; }
1050
1051   /**
1052    * Return the offset for random mnumbers.
1053    */
1054   double offset() const { return theOffset; }
1055
1056   /**
1057    * Return the density's value
1058    */
1059   double value(double x) const {
1060     return x>=lower() && x<=upper() ? 1/abs(x-pole()) : 0.;
1061   }
1062
1063   /**
1064    * Return the density's normalization
1065    */
1066   double normalization() const {
1067     return scale();
1068   }
1069
1070   /**
1071    * Generate the return value according to the implemented density,
1072    * given a flat random number on the unit interval.
1073    */
1074   double operator()(double r) const {
1075     return pole() + sign(upper()-pole())*exp(scale()*r+offset());
1076   }
1077
1078 };
1079
1080 /**
1081  * \ingroup Matchbox
1082  * \author Simon Platzer
1083  * \brief Construct the density 1/|x-z|
1084  */
1085 inline Generator<Inverse> inverse(double z,
1086                   double lower, double upper) {
1087   return Generator<Inverse>(z,lower,upper);
1088 }
1089
1090 /**
1091  * \ingroup Matchbox
1092  * \author Simon Platzer
1093  * \brief The density |(x-z)|^p
1094  */
1095 struct Power {};
1096
1097 /**
1098  * \ingroup Matchbox
1099  * \author Simon Platzer
1100  * \brief Generate x with density |(x-z)|^p
1101  */
1102 template<>
1103 class Generator<Power> {
1104
1105   /**
1106    * The position of the pole
1107    */
1108   double thePole;
1109
1110   /**
1111    * The power
1112    */
1113   double thePower;
1114
1115   /**
1116    * The lower bound
1117    */
1118   double theLower;
1119
1120   /**
1121    * The upper bound
1122    */
1123   double theUpper;
1124
1125   /**
1126    * Scale for random numbers.
1127    */
1128   double theScale;
1129
1130   /**
1131    * Offset for random mnumbers.
1132    */
1133   double theOffset;
1134
1135 public:
1136
1137   /**
1138    * Construct from pole, power and boundaries.
1139    */
1140   Generator(double z, double p,
1141         double l, double u)
1142     : thePole(z), thePower(p),
1143       theLower(l), theUpper(u),
1144       theScale(z<=l ? (pow(u-z,1.+p)-pow(l-z,1.+p))/(1.+p) : (pow(z-l,1.+p)-pow(z-u,1.+p))/(1.+p)),
1145       theOffset(z<=l ? pow(l-z,1.+p)/(1.+p) : pow(z-u,1.+p)/(1.+p)) {
1146     if ( p == -1. )
1147       throw std::logic_error("[Generator<Power>] Unit inverse. Consider using inverse().");
1148     if ( z >= l && z <= u && p < 0. )
1149       throw std::logic_error("[Generator<Power>] Pole inside sampling interval.");
1150     if ( z >= l && z <= u && p > 0. )
1151       throw std::logic_error("[Generator<Power>] Zero inside sampling interval.");
1152   }
1153
1154   /**
1155    * Return the lower bound of the density generated.
1156    */
1157   double lower() const { return theLower; }
1158
1159   /**
1160    * Return the upper bound of the density generated.
1161    */
1162   double upper() const { return theUpper; }
1163
1164   /**
1165    * Return the position of the pole.
1166    */
1167   double pole() const { return thePole; }
1168
1169   /**
1170    * Return the power.
1171    */
1172   double power() const { return thePower; }
1173
1174   /**
1175    * Return the scale for random numbers.
1176    */
1177   double scale() const { return theScale; }
1178
1179   /**
1180    * Return the offset for random mnumbers.
1181    */
1182   double offset() const { return theOffset; }
1183
1184   /**
1185    * Return the density's value
1186    */
1187   double value(double x) const {
1188     return x>=lower() && x<=upper() ? pow(abs(x-pole()),power()) : 0.;
1189   }
1190
1191   /**
1192    * Return the density's normalization
1193    */
1194   double normalization() const {
1195     return scale();
1196   }
1197
1198   /**
1199    * Generate the return value according to the implemented density,
1200    * given a flat random number on the unit interval.
1201    */
1202   double operator()(double r) const {
1203     return pole() + sign(upper()-pole())*pow((1.+power())*(scale()*r+offset()),1./(1.+power()));
1204   }
1205
1206 };
1207
1208 /**
1209  * \ingroup Matchbox
1210  * \author Simon Platzer
1211  * \brief Construct the density |(x-z)|^p
1212  */
1213 inline Generator<Power> power(double z, double p,
1214                   double lower, double upper) {
1215   return Generator<Power>(z,p,lower,upper);
1216 }
1217
1218 /**
1219  * \ingroup Matchbox
1220  * \author Simon Platzer
1221  * \brief The density 1/((x-z)^2 + abs(w z))
1222  */
1223 struct BreitWigner {};
1224
1225 /**
1226  * \ingroup Matchbox
1227  * \author Simon Platzer
1228  * \brief Generate x with density 1/((x-z)^2 + abs(w z))
1229  */
1230 template<>
1231 class Generator<BreitWigner> {
1232
1233   /**
1234    * The position of the pole
1235    */
1236   double thePole;
1237
1238   /**
1239    * The width
1240    */
1241   double theWidth;
1242
1243   /**
1244    * The lower bound
1245    */
1246   double theLower;
1247
1248   /**
1249    * The upper bound
1250    */
1251   double theUpper;
1252
1253   /**
1254    * Scale for random numbers.
1255    */
1256   double theScale;
1257
1258   /**
1259    * Offset for random mnumbers.
1260    */
1261   double theOffset;
1262
1263   /**
1264    * The square root of width times pole.
1265    */
1266   double theSqrtWidth;
1267
1268 public:
1269
1270   /**
1271    * Construct from pole, width and boundaries.
1272    */
1273   Generator(double z, double w,
1274         double l, double u)
1275     : thePole(z), theWidth(w),
1276       theLower(l), theUpper(u),
1277       theScale((atan((u-z)/sqrt(abs(w*z)))-atan((l-z)/sqrt(abs(w*z))))/sqrt(abs(w*z))),
1278       theOffset(atan((l-z)/sqrt(abs(w*z)))/sqrt(abs(w*z))),
1279       theSqrtWidth(sqrt(abs(w*z))) {
1280     if ( w == 0. )
1281       throw std::logic_error("[Generator<BreitWigner>] Zero width. Consider using power().");
1282   }
1283
1284   /**
1285    * Return the lower bound of the density generated.
1286    */
1287   double lower() const { return theLower; }
1288
1289   /**
1290    * Return the upper bound of the density generated.
1291    */
1292   double upper() const { return theUpper; }
1293
1294   /**
1295    * Return the position of the pole.
1296    */
1297   double pole() const { return thePole; }
1298
1299   /**
1300    * Return the width.
1301    */
1302   double width() const { return theWidth; }
1303
1304   /**
1305    * Return the scale for random numbers.
1306    */
1307   double scale() const { return theScale; }
1308
1309   /**
1310    * Return the offset for random mnumbers.
1311    */
1312   double offset() const { return theOffset; }
1313
1314   /**
1315    * The square root of width times pole.
1316    */
1317   double sqrtWidth() const { return theSqrtWidth; }
1318
1319   /**
1320    * Return the density's value
1321    */
1322   double value(double x) const {
1323     return
1324       x>=lower() && x<=upper() ?
1325       1./(sqr(x-pole())+abs(width()*pole())) : 0.;
1326   }
1327
1328   /**
1329    * Return the density's normalization
1330    */
1331   double normalization() const {
1332     return scale();
1333   }
1334
1335   /**
1336    * Generate the return value according to the implemented density,
1337    * given a flat random number on the unit interval.
1338    */
1339   double operator()(double r) const {
1340     double res = pole() + sqrtWidth()*tan(sqrtWidth()*(scale()*r+offset()));
1341     if ( res <= lower() ) return lower()*(1+std::numeric_limits<double>::epsilon());
1342     else if ( res >= upper() ) return upper()*(1-std::numeric_limits<double>::epsilon());
1343     else return res;
1344   }
1345
1346 };
1347
1348 /**
1349  * \ingroup Matchbox
1350  * \author Simon Platzer
1351  * \brief Construct the density 1/((x-z)^2 + abs(w z))
1352  */
1353 inline Generator<BreitWigner> breitWigner(double z, double w,
1354                       double lower, double upper) {
1355   return Generator<BreitWigner>(z,w,lower,upper);
1356 }
1357
1358 }
1359
1360 }
1361
1362 #endif // HERWIG_RandomHelpers_H