ThePEG/Handlers/ACDCSampler.h

0001 // -*- C++ -*-
0002 //
0003 // ACDCSampler.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 1999-2019 Leif Lonnblad
0005 //
0006 // ThePEG 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 ThePEG_ACDCSampler_H
0010 #define ThePEG_ACDCSampler_H
0011 // This is the declaration of the ACDCSampler class.
0012
0013 #include "ThePEG/Handlers/SamplerBase.h"
0014 #include "ThePEG/ACDC/ACDCGen.h"
0015 #include "ThePEG/Handlers/StandardEventHandler.h"
0016 #include "ThePEG/Repository/RandomGenerator.h"
0017 #include "ThePEG/Repository/UseRandom.h"
0018 #include "ThePEG/Utilities/SimplePhaseSpace.xh"
0019
0020 // #include "ACDCSampler.fh"
0021 // #include "ACDCSampler.xh"
0022
0023 namespace ThePEG {
0024
0025 /**
0026  * This class inherits from SampleBase and implements
0027  * the Auto Compensating Divide-and-Conquer phase space generator,
0028  * ACDCGen.
0029  *
0030  * @see \ref ACDCSamplerInterfaces "The interfaces"
0031  * defined for ACDCSampler.
0032  * @see ACDCGen
0033  */
0034 class ACDCSampler: public SamplerBase {
0035
0036 public:
0037
0038   /** Typedef the underlying ACDCGen class. */
0039   typedef ACDCGenerator::ACDCGen<UseRandom,tStdEHPtr> SamplerType;
0040
0041   /** @name Standard constructors and destructors. */
0042   //@{
0043   /**
0044    * The default constructor.
0045    */
0046   ACDCSampler() : theEps(100*Constants::epsilon), theMargin(1.1), theNTry(1000) {}
0047
0048   /**
0049    * The copy constructor. We don't copy theSampler.
0050    */
0051   ACDCSampler(const ACDCSampler & x)
0052     : SamplerBase(x), theSampler(),
0053       theEps(x.theEps), theMargin(x.theMargin),
0054       theNTry(x.theNTry) {}
0055
0056   /**
0057    * The destructor.
0058    */
0059   virtual ~ACDCSampler();
0060   //@}
0061
0062 public:
0063
0064   /** @name Virtual functions needed for SamplerBase */
0065   //@{
0066   /**
0067    * Initialize the the sampler, possibly doing presampling of the
0068    * phase space.
0069    */
0070   virtual void initialize();
0071
0072   /**
0073    * Generarate a new phase space point and return a weight associated
0074    * with it. This weight should preferably be 1.
0075    */
0076   virtual double generate();
0077
0078   /**
0079    * ACDCSampler is able to sample several different functions
0080    * separately. This function returns the last chosen
0081    * function.
0082    */
0083   virtual int lastBin() const;
0084
0085   /**
0086    * Reject the last chosen phase space point.
0087    */
0088   virtual void rejectLast();
0089
0090   /**
0091    * Return the total integrated cross section determined from the
0092    * Monte Carlo sampling so far.
0093    */
0094   virtual CrossSection integratedXSec() const;
0095
0096   /**
0097    * Return the error on the total integrated cross section determined
0098    * from the Monte Carlo sampling so far.
0099    */
0100   virtual CrossSection integratedXSecErr() const;
0101
0102   /**
0103    * Return the sum of the weights returned by generate() so far (of
0104    * the events that were not rejeted).
0105    */
0106   virtual double sumWeights() const;
0107
0108   /**
0109    * Return the sum of the weights squared returned by generate() so far (of
0110    * the events that were not rejeted).
0111    */
0112   virtual double sumWeights2() const;
0113   //@}
0114
0115 public:
0116
0117   /** @name Functions used by the persistent I/O system. */
0118   //@{
0119   /**
0120    * Function used to write out object persistently.
0121    * @param os the persistent output stream written to.
0122    */
0123   void persistentOutput(PersistentOStream & os) const;
0124
0125   /**
0126    * Function used to read in object persistently.
0127    * @param is the persistent input stream read from.
0128    * @param version the version number of the object when written.
0129    */
0130   void persistentInput(PersistentIStream & is, int version);
0131   //@}
0132
0133   /**
0134    * The standard Init function used to initialize the interfaces.
0135    * Called exactly once for each class by the class description system
0136    * before the main function starts or
0137    * when this class is dynamically loaded.
0138    */
0139   static void Init();
0140
0141 protected:
0142
0143   /** @name Clone Methods. */
0144   //@{
0145   /**
0146    * Make a simple clone of this object.
0147    * @return a pointer to the new object.
0148    */
0149   virtual IBPtr clone() const;
0150
0151   /** Make a clone of this object, possibly modifying the cloned object
0152    * to make it sane.
0153    * @return a pointer to the new object.
0154    */
0155   virtual IBPtr fullclone() const;
0156
0157   //@}
0158
0159 protected:
0160
0161   /** @name Standard Interfaced functions. */
0162   //@{
0163   /**
0164    * Initialize this object. Called in the run phase just before
0165    * a run begins.
0166    */
0167   virtual void doinitrun();
0168
0169   /**
0170    * Finalize this object. Called in the run phase just after a
0171    * run has ended. Used eg. to write out statistics.
0172    */
0173   virtual void dofinish();
0174   //@}
0175
0176 private:
0177
0178   /**
0179    * The actual sampler object.
0180    */
0181   SamplerType theSampler;
0182
0183   /**
0184    * The smallest possible division allowed.
0185    */
0186   double theEps;
0187
0188   /**
0189    * The factor controlling the loss of efficiency when compensating.
0190    */
0191   double theMargin;
0192
0193   /**
0194    * The number of points to use to find initial average.
0195    */
0196   int theNTry;
0197
0198 protected:
0199
0200   /** @cond EXCEPTIONCLASSES */
0201   /** Exception class used by ACDCSampler if the undelying ACDCGen was
0202       still in a compensating mode when the run was finished */
0203   struct ACDCStillCompensating: public Exception {};
0204
0205   /** Exception class used by ACDCSampler if a StandardEventHandler
0206       was not able to produce a non-zero cross section. */
0207   struct EventInitNoXSec: public InitException {};
0208
0209   /** Exception class used if ACDCSampler was not able to produce a
0210       phase space point within the maximum allowed number of
0211       attempts. */
0212   struct EventLoopException: public Exception {};
0213   /** @endcond */
0214
0215 private:
0216
0217   /**
0218    * Describe a concrete class with persistent data.
0219    */
0220   static ClassDescription<ACDCSampler> initACDCSampler;
0221
0222   /**
0223    *  Private and non-existent assignment operator.
0224    */
0225   ACDCSampler & operator=(const ACDCSampler &) = delete;
0226
0227 };
0228
0229 }
0230
0231 namespace ThePEG {
0232
0233 /** @cond TRAITSPECIALIZATIONS */
0234
0235 /**
0236  * The following template specialization informs ThePEG about the
0237  * base class of ACDCSampler.
0238  */
0239 template <>
0240 struct BaseClassTrait<ACDCSampler,1>: public ClassTraitsType {
0241   /** Typedef of the first base class of ACDCSampler. */
0242   typedef SamplerBase NthBase;
0243 };
0244
0245 /**
0246  * The following template specialization informs ThePEG about the
0247  * name of this class and the shared object where it is defined.
0248  */
0249 template <>
0250 struct ClassTraits<ACDCSampler>: public ClassTraitsBase<ACDCSampler> {
0251   /**
0252    * Return the class name.
0253    */
0254   static string className() { return "ThePEG::ACDCSampler"; }
0255   /** Return the name of the shared library to be loaded to get
0256    * access to this class and every other class it uses
0257    * (except the base class).
0258    */
0259   static string library() { return "ACDCSampler.so"; }
0260
0261 };
0262
0263 /** @endcond */
0264
0265 }
0266
0267 namespace ACDCGenerator {
0268
0269 /** @cond TRAITSPECIALIZATIONS */
0270
0271 /** Specialized Traits class to define the interface to the
0272  * StandardEventHandler object to be sampled by ACDCGen.
0273  */
0274 template <>
0275 struct ACDCFncTraits<ThePEG::tStdEHPtr>: public ACDCTraitsType {
0276   /** Convenient typdef. */
0277   typedef ThePEG::tStdEHPtr tStdEHPtr;
0278   /**
0279    * Call a function to be sampled by ACDCGen.
0280    * @return <code>(*f)(x)</code>.
0281    */
0282   static inline double value(const tStdEHPtr & eh, const DVector & x) {
0283     using namespace ThePEG::Units;
0284     try {
0285       return eh->dSigDR(x)/nanobarn;
0286     }
0287     catch ( ThePEG::ImpossibleKinematics & v ) {
0288       breakThePEG();
0289     }
0290     catch ( std::exception & e ) {
0291       breakThePEG();
0292     }
0293     catch ( ... ) {
0294       breakThePEG();
0295     }
0296     return 0.0;
0297   }
0298
0299 };
0300
0301 /** Specialized Traits class to inform ACDCGen how to use the
0302     static UseRandom class. */
0303 template <>
0304 struct ACDCRandomTraits<ThePEG::UseRandom>: public ACDCTraitsType {
0305   /** Convenient typedef. */
0306   typedef ThePEG::UseRandom UseRandom;
0307
0308   /**
0309    * Return a flat random number in the interval ]0,1[.
0310    */
0311   static inline double rnd(UseRandom *) { return UseRandom::rnd(); }
0312
0313   /**
0314    * Return a flat random number in the interval ]\a xl,\a xu[.
0315    */
0316   static inline double rnd(UseRandom * r, double xl, double xu) {
0317     return xl + (xu - xl)*rnd(r);
0318   }
0319
0320   /**
0321    * Generate a set of random numbers.
0322    * @param r the random generator.
0323    * @param l an input iterator giving the lower limit of the interval
0324    * of the first requested random number.
0325    * @param lend an input iterator marking the end of the range of
0326    * requested random numbers.
0327    * @param u an input iterator giving the upper limit of the interval
0328    * of the first requested random number.
0329    * @param res the ouput iterator used to output the random numbers.
0330    */
0331   template <typename InputIterator, typename OutputIterator>
0332   static inline void rnd(UseRandom * r,
0333              InputIterator l, InputIterator lend,
0334              InputIterator u, OutputIterator res) {
0335     for ( ; l != lend; ++l ) *res++ = *l + (*u++ - *l)*rnd(r);
0336   }
0337
0338   /**
0339    * Generate \a D random numbers. The numbers are put into the
0340    * OutputIterator \a res.
0341    */
0342   template <typename OutputIterator>
0343   static inline void rnd(UseRandom * r, int D, OutputIterator res) {
0344     for ( int d = 0; d < D; ++d ) *res++ = rnd(r);
0345   }
0346
0347   /**
0348    * Return true with probability \a x.
0349    */
0350   static inline bool rndBool(UseRandom, double x) {
0351     return UseRandom::rndbool(x);
0352   }
0353
0354   /**
0355    * Return true with probability \a x(\a x + \a y).
0356    */
0357   static inline bool rndBool(UseRandom *, double x, double y) {
0358     return UseRandom::rndbool(x, y);
0359   }
0360
0361   /**
0362    * Return a random integer in the interval [0,\a x[.
0363    */
0364   static inline long rndInt(UseRandom *, long x) {
0365     return UseRandom::irnd(x);
0366   }
0367
0368 };
0369
0370 /** @endcond */
0371
0372 }
0373
0374 #endif /* ThePEG_ACDCSampler_H */