DIRE/Shower/Shower.H

0001 #ifndef DIRE__Shower__Shower_H
0002 #define DIRE__Shower__Shower_H
0003
0004 #include "DIRE/Tools/Parton.H"
0005 #include "DIRE/Shower/Kernel.H"
0006 #include "DIRE/Tools/Amplitude.H"
0007 #include "ATOOLS/Phys/Weights.H"
0008
0009 namespace ATOOLS {
0010   class QCD_Variation_Params;
0011 }
0012 namespace MODEL {
0013   class Model_Base;
0014   class Running_AlphaS;
0015 }
0016 namespace PDF {
0017   class PDF_Base;
0018   class ISR_Handler;
0019 }
0020
0021 namespace DIRE {
0022
0023   class Cluster_Definitions;
0024
0025   class Shower {
0026   public:
0027
0028     struct Reweight_Args {
0029       Splitting *m_s;
0030       int m_acc;
0031       Reweight_Args(Splitting *const s,const int acc):
0032     m_s(s), m_acc(acc) {}
0033     };// end of struct Reweight_Args
0034
0035     struct JetVeto_Args {
0036       ATOOLS::Cluster_Amplitude *p_ampl;
0037       double m_jcv;
0038       int m_acc, m_nlo;
0039       std::vector<int> m_skip;
0040       JetVeto_Args(ATOOLS::Cluster_Amplitude *const ampl,
0041            const double &jcv,const size_t &n):
0042     p_ampl(ampl), m_jcv(jcv), m_acc(0), m_skip(n,0) {}
0043     };// end of struct JetVeto_Args
0044
0045     typedef std::map<ATOOLS::Flavour,Kernel_Vector> SKernel_Map;
0046
0047     typedef std::map<ATOOLS::Flavour,Kernel*> EKernel_Map;
0048     typedef std::map<ATOOLS::Flavour,EKernel_Map> SEKernel_Map;
0049     typedef std::map<int,SEKernel_Map> Kernel_Map;
0050
0051   private:
0052
0053     // set up a C array wrapper used to hold a table of cumulative integrals in
0054     // GeneratePoint
0055     // rows: kernels allowed by the splitter
0056     // cols: possible spectators
0057     // NOTE: by using (encapsulated) C arrays we can avoid using vectors of
0058     // vectors in the hot-spot function GeneratePoint
0059     struct CumulativeIntegralTable {
0060       CumulativeIntegralTable() : nrows {0}, ncols {0} {};
0061       CumulativeIntegralTable(int _nrows, int _ncols) :
0062         nrows{_nrows}, ncols {_ncols}
0063       {
0064         sumsizes = new int[nrows];
0065         spects = new size_t[nrows*ncols];
0066         sums = new double[nrows*ncols];
0067       }
0068       CumulativeIntegralTable& operator=(CumulativeIntegralTable&& other)
0069       {
0070         if (this == &other)
0071           return *this;
0072         if (nrows > 0) {
0073           delete[] sumsizes;
0074           delete[] spects;
0075           delete[] sums;
0076         }
0077         nrows = other.nrows;
0078         ncols = other.ncols;
0079         sumsizes = other.sumsizes;
0080         spects = other.spects;
0081         sums = other.sums;
0082         other.nrows = 0;
0083         return *this;
0084       }
0085       ~CumulativeIntegralTable() {
0086         if (nrows > 0) {
0087           delete[] sumsizes;
0088           delete[] spects;
0089           delete[] sums;
0090         }
0091       }
0092       void Clear(int row) { sumsizes[row] = 0; }
0093       int Size(int row) { return sumsizes[row]; }
0094       double Sum(int row, int col) { return sums[row*ncols + col]; }
0095       size_t Spect(int row, int col) { return spects[row*ncols + col]; }
0096       double LastSum(int row) { return sums[row*ncols + sumsizes[row] - 1]; }
0097       double LastSpect(int row) { return spects[row*ncols + sumsizes[row] - 1]; }
0098       void AppendSumAndSpect(int row, double sum, size_t spect) {
0099         const auto idx = row*ncols + sumsizes[row];
0100         sums[idx] = sum;
0101         spects[idx] = spect;
0102         ++(sumsizes[row]);
0103       }
0104       int nrows, ncols;
0105       int* sumsizes;
0106       size_t* spects;
0107       double* sums;
0108     };
0109
0110     CumulativeIntegralTable m_sums;
0111
0112     MODEL::Model_Base *p_model;
0113     PDF::PDF_Base     *p_pdf[2];
0114
0115     MODEL::Running_AlphaS *p_as;
0116
0117     Cluster_Definitions *p_cluster;
0118
0119     ATOOLS::Weights_Map m_weightsmap;
0120     bool m_reweight;
0121
0122     Kernel_Vector m_cks;
0123     SKernel_Map   m_sks;
0124     Kernel_Map    m_kmap;
0125
0126     double m_tmin[2], m_rewtmin, m_cplfac[2], m_rsf, m_fsf;
0127     double m_weight, m_oef, m_pdfmin[2], m_rcf, m_tcef;
0128
0129     int m_kin, m_kfac, m_cpl, m_mec;
0130
0131     unsigned int m_maxem, m_maxpart;
0132
0133     void AddKernel(Kernel *const k);
0134
0135     void AddWeight(const Amplitude &a,const double &t);
0136
0137     Splitting GeneratePoint(Parton &p,const double &t,
0138                 const int &cm,const unsigned int &nem);
0139     Splitting GeneratePoint(const Amplitude &a,const double &t,
0140                 const unsigned int &nem);
0141
0142     void Reweight(ATOOLS::QCD_Variation_Params* params,
0143                   size_t varindex,
0144                   const Reweight_Args& a);
0145
0146   public:
0147
0148     Shower();
0149
0150     ~Shower();
0151
0152     bool Init(MODEL::Model_Base *const,
0153           PDF::ISR_Handler *const);
0154
0155     void SetMS(const ATOOLS::Mass_Selector *const ms);
0156
0157     int Evolve(Amplitude&, unsigned int& nem);
0158
0159     double GetXPDF(const double &x,const double &Q2,
0160            const ATOOLS::Flavour &fl,const int b) const;
0161
0162     Kernel *GetKernel(const Splitting &s,const int mode) const;
0163
0164     int RemnantTest(Parton *const c,const ATOOLS::Vec4D &p);
0165
0166     inline MODEL::Model_Base *Model() const { return p_model; }
0167
0168     inline MODEL::Running_AlphaS *const &AlphaS() const { return p_as; }
0169
0170     inline double TMin(const int i) const { return m_tmin[i]; }
0171
0172     inline double CplFac(const int i) const { return m_cplfac[i]; }
0173
0174     inline double PDFMin(const int i) const { return m_pdfmin[i]; }
0175
0176     inline double GetWeight() const { return m_weight; }
0177     inline const ATOOLS::Weights_Map& GetWeightsMap() const
0178     {
0179       return m_weightsmap;
0180     }
0181
0182     inline int KFactorScheme() const  { return m_kfac; }
0183     inline int CouplingScheme() const { return m_cpl;  }
0184
0185     inline int KinematicsScheme() const { return m_kin; }
0186
0187     inline int MECorrection() const { return m_mec; }
0188
0189     inline double MuR2Factor() const { return m_rsf; }
0190
0191   };//end of class Shower
0192
0193 }// end of namespace DIRE
0194
0195 #endif