Pythia8Interface/plugins/ReweightUserHooks.h

0001 #include "Pythia8/Pythia.h"
0002 #include "TF1.h"
0003
0004 class PtHatReweightUserHook : public Pythia8::UserHooks {
0005 public:
0006   PtHatReweightUserHook(double _pt = 15, double _power = 4.5) : pt(_pt), power(_power) {}
0007   ~PtHatReweightUserHook() override {}
0008
0009   bool canBiasSelection() override { return true; }
0010
0011   double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr,
0012                          const Pythia8::PhaseSpace* phaseSpacePtr,
0013                          bool inEvent) override {
0014     //the variable selBias of the base class should be used;
0015     if ((sigmaProcessPtr->nFinal() == 2)) {
0016       selBias = pow(phaseSpacePtr->pTHat() / pt, power);
0017       return selBias;
0018     }
0019     selBias = 1.;
0020     return selBias;
0021   }
0022
0023 private:
0024   double pt, power;
0025 };
0026
0027 class PtHatEmpReweightUserHook : public Pythia8::UserHooks {
0028 public:
0029   PtHatEmpReweightUserHook(const std::string& tuneName = "") {
0030     if (tuneName == "CP5" || tuneName == "CP5Run3")
0031       p = {7377.94700788, 8.38168461349, -4.70983112392, -0.0310148108446, -0.028798537937, 925.335472326};
0032     //Default reweighting - works good for tune CUEPT8M1
0033     else
0034       p = {5.3571961909810e+13,
0035            1.0907678218282e+01,
0036            -2.5898069229451e+00,
0037            -5.1575514014931e-01,
0038            5.5951279807561e-02,
0039            3.5e+02};
0040     const double ecms = (tuneName == "CP5Run3" ? 13600. : 13000.);
0041     sigma = [this, ecms](double x) -> double {
0042       return (p[0] * pow(x, p[2] + p[3] * log(0.01 * x) + p[4] * pow(log(0.01 * x), 2)) *
0043               pow(1 - 2 * x / (ecms + p[5]), p[1])) *
0044              x;
0045     };
0046   }
0047   ~PtHatEmpReweightUserHook() override {}
0048
0049   bool canBiasSelection() override { return true; }
0050
0051   double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr,
0052                          const Pythia8::PhaseSpace* phaseSpacePtr,
0053                          bool inEvent) override {
0054     //the variable selBias of the base class should be used;
0055     if ((sigmaProcessPtr->nFinal() == 2)) {
0056       selBias = 1.0 / sigma(phaseSpacePtr->pTHat());
0057       return selBias;
0058     }
0059     selBias = 1.;
0060     return selBias;
0061   }
0062
0063 private:
0064   std::vector<double> p;
0065   std::function<double(double)> sigma;
0066 };
0067
0068 class RapReweightUserHook : public Pythia8::UserHooks {
0069 public:
0070   RapReweightUserHook(const std::string& _yLabsigma_func,
0071                       double _yLab_power,
0072                       const std::string& _yCMsigma_func,
0073                       double _yCM_power,
0074                       double _pTHatMin,
0075                       double _pTHatMax)
0076       : yLabsigma_func(_yLabsigma_func),
0077         yCMsigma_func(_yCMsigma_func),
0078         yLab_power(_yLab_power),
0079         yCM_power(_yCM_power),
0080         pTHatMin(_pTHatMin),
0081         pTHatMax(_pTHatMax) {
0082     // empirical parametrizations defined in configuration file
0083     yLabsigma = TF1("yLabsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
0084     yCMsigma = TF1("yCMsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
0085   }
0086   ~RapReweightUserHook() override {}
0087
0088   bool canBiasSelection() override { return true; }
0089
0090   double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr,
0091                          const Pythia8::PhaseSpace* phaseSpacePtr,
0092                          bool inEvent) override {
0093     //the variable selBias of the base class should be used;
0094     if ((sigmaProcessPtr->nFinal() == 2)) {
0095       double x1 = phaseSpacePtr->x1();
0096       double x2 = phaseSpacePtr->x2();
0097       double yLab = 0.5 * log(x1 / x2);
0098       double yCM = 0.5 * log(phaseSpacePtr->tHat() / phaseSpacePtr->uHat());
0099       double pTHat = phaseSpacePtr->pTHat();
0100       double sigmaLab = yLabsigma.Eval(pTHat);
0101       double sigmaCM = yCMsigma.Eval(pTHat);
0102       // empirical reweighting function
0103       selBias = exp(pow(fabs(yLab), yLab_power) / (2 * sigmaLab * sigmaLab) +
0104                     pow(fabs(yCM), yCM_power) / (2 * sigmaCM * sigmaCM));
0105       return selBias;
0106     }
0107     selBias = 1.;
0108     return selBias;
0109   }
0110
0111 private:
0112   std::string yLabsigma_func, yCMsigma_func;
0113   double yLab_power, yCM_power, pTHatMin, pTHatMax;
0114   TF1 yLabsigma, yCMsigma;
0115 };
0116
0117 class PtHatRapReweightUserHook : public Pythia8::UserHooks {
0118 public:
0119   PtHatRapReweightUserHook(const std::string& _yLabsigma_func,
0120                            double _yLab_power,
0121                            const std::string& _yCMsigma_func,
0122                            double _yCM_power,
0123                            double _pTHatMin,
0124                            double _pTHatMax,
0125                            double _pt = 15,
0126                            double _power = 4.5)
0127       : yLabsigma_func(_yLabsigma_func),
0128         yCMsigma_func(_yCMsigma_func),
0129         yLab_power(_yLab_power),
0130         yCM_power(_yCM_power),
0131         pTHatMin(_pTHatMin),
0132         pTHatMax(_pTHatMax),
0133         pt(_pt),
0134         power(_power) {
0135     // empirical parametrizations defined in configuration file
0136     yLabsigma = TF1("yLabsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
0137     yCMsigma = TF1("yCMsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
0138   }
0139   ~PtHatRapReweightUserHook() override {}
0140
0141   bool canBiasSelection() override { return true; }
0142
0143   double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr,
0144                          const Pythia8::PhaseSpace* phaseSpacePtr,
0145                          bool inEvent) override {
0146     //the variable selBias of the base class should be used;
0147     if ((sigmaProcessPtr->nFinal() == 2)) {
0148       double x1 = phaseSpacePtr->x1();
0149       double x2 = phaseSpacePtr->x2();
0150       double yLab = 0.5 * log(x1 / x2);
0151       double yCM = 0.5 * log(phaseSpacePtr->tHat() / phaseSpacePtr->uHat());
0152       double pTHat = phaseSpacePtr->pTHat();
0153       double sigmaLab = yLabsigma.Eval(pTHat);
0154       double sigmaCM = yCMsigma.Eval(pTHat);
0155       // empirical reweighting function
0156       selBias = pow(pTHat / pt, power) * exp(pow(fabs(yLab), yLab_power) / (2 * sigmaLab * sigmaLab) +
0157                                              pow(fabs(yCM), yCM_power) / (2 * sigmaCM * sigmaCM));
0158       return selBias;
0159     }
0160     selBias = 1.;
0161     return selBias;
0162   }
0163
0164 private:
0165   std::string yLabsigma_func, yCMsigma_func;
0166   double yLab_power, yCM_power, pTHatMin, pTHatMax, pt, power;
0167   TF1 yLabsigma, yCMsigma;
0168 };