ThePEG/ACDC/ACDCGenCell.icc

0001 // -*- C++ -*-
0002 //
0003 // ACDCGenCell.icc 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
0010 namespace ACDCGenerator {
0011
0012 inline ACDCGenCell::ACDCGenCell(double newG)
0013   : theG(newG), theV(1.0), theUpper(0), theLower(0),
0014     theDivision(-1.0), theSplitDimension(-1) {}
0015
0016 inline ACDCGenCell::
0017 ACDCGenCell(double newG, double newV)
0018   : theG(newG), theV(newV), theUpper(0), theLower(0),
0019     theDivision(-1.0), theSplitDimension(-1) {}
0020
0021 inline ACDCGenCell::~ACDCGenCell() {
0022   if ( isSplit() ) {
0023     delete theUpper;
0024     delete theLower;
0025   }
0026 }
0027
0028 template <typename RndType>
0029 inline ACDCGenCell * ACDCGenCell::
0030 generate(DVector & lo, DVector & up, RndType * rnd) {
0031   if ( isSplit() ) {
0032     if ( ACDCRandomTraits<RndType>::
0033      rndBool(rnd, upper()->maxInt(), lower()->maxInt()) ) {
0034       lo[dim()] = div();
0035       return upper()->generate(lo, up, rnd);
0036     } else {
0037       up[dim()] = div();
0038       return lower()->generate(lo, up, rnd);
0039     }
0040   }
0041   return this;
0042 }
0043
0044 inline ACDCGenCell * ACDCGenCell::
0045 generate(DVector & lo, DVector & up, DVector & rndv) {
0046   if ( isSplit() ) {
0047     double r = lower()->maxInt()/maxInt();
0048     if ( rndv[dim()] > r ) {
0049       lo[dim()] = div();
0050       rndv[dim()] = (rndv[dim()] - r)/(1.0 - r);
0051       return upper()->generate(lo, up, rndv);
0052     } else {
0053       up[dim()] = div();
0054       rndv[dim()] = rndv[dim()]/r;
0055       return lower()->generate(lo, up, rndv);
0056     }
0057   }
0058   return this;
0059 }
0060
0061 inline ACDCGenCell *  ACDCGenCell::
0062 getCell(DVector & lo, const DVector & x, DVector & up) {
0063   if ( isSplit() ) {
0064     if ( x[dim()] > div() ) {
0065       lo[dim()] = div();
0066       return upper()->getCell(lo, x, up);
0067     } else {
0068       up[dim()] = div();
0069       return lower()->getCell(lo, x, up);
0070     }
0071   } else
0072     return this;
0073 }
0074
0075 inline void ACDCGenCell::
0076 splitme(double lo, double newDiv, double up, DimType newDim) {
0077   theSplitDimension = newDim;
0078   theDivision = newDiv;
0079   double frac = (up - div())/(up - lo);
0080   theUpper = new ACDCGenCell(g(), v()*frac);
0081   frac = (div() - lo)/(up - lo);
0082   theLower = new ACDCGenCell(g(), v()*frac);
0083 }
0084
0085 inline bool ACDCGenCell::isSplit() const {
0086   return upper();
0087 }
0088
0089 inline double ACDCGenCell::doMaxInt(double scale) {
0090   if ( isSplit() ) theG = (upper()->doMaxInt() + lower()->doMaxInt())/v();
0091   else theG *= scale;
0092   return maxInt();
0093 }
0094
0095 inline void ACDCGenCell::smooth(double frac) {
0096   if ( !isSplit() ) return;
0097   upper()->smooth(frac);
0098   lower()->smooth(frac);
0099   if ( upper()->maxInt() < lower()->maxInt()*frac &&
0100        upper()->maxInt() > 0.0 )
0101     upper()->doMaxInt(lower()->maxInt()*frac/upper()->maxInt());
0102   else if ( lower()->maxInt() < upper()->maxInt()*frac &&
0103         lower()->maxInt() > 0.0 )
0104     lower()->doMaxInt(upper()->maxInt()*frac/lower()->maxInt());
0105   doMaxInt();
0106 }
0107
0108 inline double ACDCGenCell::maxInt() const {
0109   return g()*v();
0110 }
0111
0112 inline void ACDCGenCell::g(double newG) {
0113   theG = newG;
0114 }
0115
0116 inline int ACDCGenCell::nBins() const {
0117   return isSplit()? upper()->nBins() + lower()->nBins(): 1;
0118 }
0119
0120 inline int ACDCGenCell::depth() const {
0121   return isSplit()? std::max(upper()->depth(), lower()->depth()) + 1: 1;
0122 }
0123
0124 inline double ACDCGenCell::g() const {
0125   return theG;
0126 }
0127
0128 inline double ACDCGenCell::v() const {
0129   return theV;
0130 }
0131
0132 inline double ACDCGenCell::div() const  {
0133   return theDivision;
0134 }
0135
0136 inline DimType ACDCGenCell::dim() const {
0137   return theSplitDimension;
0138 }
0139
0140 inline ACDCGenCell * ACDCGenCell::upper() const {
0141   return theUpper;
0142 }
0143
0144 inline ACDCGenCell * ACDCGenCell::lower() const {
0145   return theLower;
0146 }
0147
0148 inline void ACDCGenCell::
0149 extract(DVector & lo, DVector & up, vector<ACDCGenCellInfo> & out) const {
0150
0151   // First add this cell which gets index isave.
0152   ACDCGenCellInfo::Index isave = out.size();
0153   out.push_back(ACDCGenCellInfo());
0154   out.back().g = g();
0155   out.back().v = v();
0156   out.back().up = up;
0157   out.back().lo = lo;
0158   out.back().iup = 0;
0159   out.back().ilo = 0;
0160
0161   if ( isSplit() ) {
0162     // If split add the upper cell.
0163     out[isave].iup = out.size();
0164     double save = lo[dim()];
0165     lo[dim()] = div();
0166     upper()->extract(lo, up, out);
0167     lo[dim()] = save;
0168
0169     // Then add the lower cell.
0170     out[isave].ilo = out.size();
0171     save = up[dim()];
0172     up[dim()] = div();
0173     lower()->extract(lo, up, out);
0174     up[dim()] = save;
0175   }
0176 }
0177
0178 template <typename OStream>
0179 inline OStream & operator<<(OStream & os, const ACDCGenCell & c) {
0180   os << c.dim() << c.div() << c.g() << c.v();
0181   if ( c.dim() < 0 ) return os;
0182   return os << *c.upper() << *c.lower();
0183 }
0184
0185 template <typename IStream>
0186 inline IStream & operator>>(IStream & is, ACDCGenCell & c) {
0187   is >> c.theSplitDimension >> c.theDivision >> c.theG >> c.theV;
0188   if ( c.dim() < 0 ) return is;
0189   c.theUpper = new ACDCGenCell(0.0);
0190   c.theLower = new ACDCGenCell(0.0);
0191   return is >> *c.theUpper >> *c.theLower;
0192 }
0193
0194 inline long ACDCGenCell::getIndex(const ACDCGenCell * c) const {
0195   long indx = -1;
0196   return getIndex(c, indx);
0197 }
0198
0199
0200 inline long ACDCGenCell::getIndex(const ACDCGenCell * c, long & indx) const {
0201   ++indx;
0202   if ( c == this ) return indx;
0203   if ( isSplit() ) {
0204     long i = upper()->getIndex(c, indx);
0205     if ( i >= 0 ) return i;
0206     return lower()->getIndex(c, indx);
0207   }
0208   return -1;
0209 }
0210
0211 inline ACDCGenCell * ACDCGenCell::getCell(long i) {
0212   long indx = -1;
0213   return getCell(i, indx);
0214 }
0215
0216 inline ACDCGenCell * ACDCGenCell::getCell(long i, long & indx) {
0217   ++indx;
0218   if ( i == indx ) return this;
0219   if ( isSplit() ) {
0220     ACDCGenCell * tmp = upper()->getCell(i, indx);
0221     if ( tmp ) return tmp;
0222     return lower()->getCell(i, indx);
0223   }
0224   return 0;
0225 }
0226
0227 }