Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-08 08:16:15

0001 // This is file contains the code to generate a dark sector shower in a strongly coupled, quasi-conformal hidden valley, often
0002 // referred to as "soft unclustered energy patterns (SUEP)" or "softbomb" events. The shower is generated in its rest frame and
0003 // for a realistic simulation this class needs to be interfaced with an event generator such as madgraph, pythia or herwig.
0004 //
0005 // The algorithm relies on arXiv:1305.5226. See arXiv:1612.00850 for a description of the model.
0006 // Please cite both papers when using this code.
0007 //
0008 //
0009 // Applicability:
0010 // The ratio of the parameters m and T (m/T) should be an O(1) number. For m/T>>1 and m/T<<1 the theoretical description of the shower is not valid.
0011 // The mass of the scalar which initiates the shower should be much larger than the mass of the mesons, in other words M>>m,T, by at least an order of magnitude.
0012 //
0013 // Written by Simon Knapen on 12/22/2019, following 1205.5226
0014 // Adapted by Carlos Erice for CMSSW
0015 #include "GeneratorInterface/Pythia8Interface/interface/SuepShower.h"
0016 
0017 // constructor
0018 SuepShower::SuepShower(double mass, double temperature, Pythia8::Rndm* rndmPtr) {
0019   edm::LogInfo("Pythia8Interface") << "Creating a SuepShower module";
0020   // Parameter setup
0021   darkmeson_mass_ = mass;
0022   fRndmPtr_ = rndmPtr;
0023   // No need to save temperature, as everything depends on it through m/T
0024   mass_over_T_ = darkmeson_mass_ / temperature;
0025 
0026   // 128 bit numerical precision, i.e. take machine precision
0027   tolerance_ = boost::math::tools::eps_tolerance<double>(128);
0028 
0029   // Median momentum
0030   p_m_ = sqrt(2 / (mass_over_T_ * mass_over_T_) * (1 + sqrt(1 + mass_over_T_ * mass_over_T_)));
0031   double pmax = sqrt(2 + 2 * sqrt(1 + mass_over_T_ * mass_over_T_)) / mass_over_T_;
0032   // Find the two cases where f(p)/f(pmax) = 1/e, given MB distribution, one appears at each side of pmax always
0033   p_plus_ = (boost::math::tools::bisect(
0034                  boost::bind(&SuepShower::logTestFunction, this, boost::placeholders::_1), pmax, 50.0, tolerance_))
0035                 .first;  // first root
0036   p_minus_ = (boost::math::tools::bisect(
0037                   boost::bind(&SuepShower::logTestFunction, this, boost::placeholders::_1), 0.0, pmax, tolerance_))
0038                  .first;  // second root
0039   // Define the auxiliar quantities for the random generation of momenta
0040   lambda_plus_ = -fMaxwellBoltzmann(p_plus_) / fMaxwellBoltzmannPrime(p_plus_);
0041   lambda_minus_ = fMaxwellBoltzmann(p_minus_) / fMaxwellBoltzmannPrime(p_minus_);
0042   q_plus_ = lambda_plus_ / (p_plus_ - p_minus_);
0043   q_minus_ = lambda_minus_ / (p_plus_ - p_minus_);
0044   q_m_ = 1 - (q_plus_ + q_minus_);
0045 }
0046 
0047 SuepShower::~SuepShower() {}
0048 
0049 // Generate a shower in the rest frame of the mediator
0050 std::vector<Pythia8::Vec4> SuepShower::generateShower(double energy) {
0051   mediator_energy_ = energy;
0052   std::vector<Pythia8::Vec4> shower;
0053   double shower_energy = 0.0;
0054 
0055   // Fill up shower record
0056   while (shower_energy < mediator_energy_) {
0057     shower.push_back(generateFourVector());
0058     shower_energy += (shower.back()).e();
0059   }
0060 
0061   // Reballance momenta to ensure conservation
0062   // Correction starts at 0
0063   Pythia8::Vec4 correction = Pythia8::Vec4(0., 0., 0., 0.);
0064   for (const auto& daughter : shower) {
0065     correction = correction + daughter;
0066   }
0067   correction = correction / shower.size();
0068   // We only want to correct momenta first
0069   correction.e(0);
0070   for (auto& daughter : shower) {
0071     daughter = daughter - correction;
0072   }
0073 
0074   // With momentum conserved, balance energy. scale is the multiplicative factor needed such that sum_daughters((scale*p)^2+m^2) = E_parent, i.e. energy is conserved
0075   double scale;
0076   scale =
0077       (boost::math::tools::bisect(
0078            boost::bind(&SuepShower::reballanceFunction, this, boost::placeholders::_1, shower), 0.0, 2.0, tolerance_))
0079           .first;
0080 
0081   for (auto& daughter : shower) {
0082     daughter.px(daughter.px() * scale);
0083     daughter.py(daughter.py() * scale);
0084     daughter.pz(daughter.pz() * scale);
0085     // Force everything to be on-shell
0086     daughter.e(sqrt(daughter.pAbs2() + darkmeson_mass_ * darkmeson_mass_));
0087   }
0088   return shower;
0089 }
0090 
0091 //////// Private Methods ////////
0092 // Maxwell-boltzman distribution, slightly massaged
0093 const double SuepShower::fMaxwellBoltzmann(double p) {
0094   return p * p * exp(-mass_over_T_ * p * p / (1 + sqrt(1 + p * p)));
0095 }
0096 
0097 // Derivative of maxwell-boltzmann
0098 const double SuepShower::fMaxwellBoltzmannPrime(double p) {
0099   return exp(-mass_over_T_ * p * p / (1 + sqrt(1 + p * p))) * p * (2 - mass_over_T_ * p * p / sqrt(1 + p * p));
0100 }
0101 
0102 // Test function to be solved for p_plus_ and p_minus_
0103 const double SuepShower::logTestFunction(double p) { return log(fMaxwellBoltzmann(p) / fMaxwellBoltzmann(p_m_)) + 1.0; }
0104 
0105 // Generate one random 4 std::vector from the thermal distribution
0106 const Pythia8::Vec4 SuepShower::generateFourVector() {
0107   double en, phi, theta, momentum;  //kinematic variables of the 4 std::vector
0108 
0109   // First do momentum, following naming of arxiv:1305.5226
0110   double U, V, X, Y, E;
0111   int i = 0;
0112   while (i < 100) {
0113     // Very rarely (<1e-5 events) takes more than one loop to converge, set to 100 for safety
0114     U = fRndmPtr_->flat();
0115     V = fRndmPtr_->flat();
0116 
0117     if (U < q_m_) {
0118       Y = U / q_m_;
0119       X = (1 - Y) * (p_minus_ + lambda_minus_) + Y * (p_plus_ - lambda_plus_);
0120       if (V < fMaxwellBoltzmann(X) / fMaxwellBoltzmann(p_m_) && X > 0) {
0121         break;
0122       }
0123     } else {
0124       if (U < q_m_ + q_plus_) {
0125         E = -log((U - q_m_) / q_plus_);
0126         X = p_plus_ - lambda_plus_ * (1 - E);
0127         if (V < exp(E) * fMaxwellBoltzmann(X) / fMaxwellBoltzmann(p_m_) && X > 0) {
0128           break;
0129         }
0130       } else {
0131         E = -log((U - (q_m_ + q_plus_)) / q_minus_);
0132         X = p_minus_ + lambda_minus_ * (1 - E);
0133         if (V < exp(E) * fMaxwellBoltzmann(X) / fMaxwellBoltzmann(p_m_) && X > 0) {
0134           break;
0135         }
0136       }
0137     }
0138   }
0139   // X is the dimensionless momentum, p/m
0140   momentum = X * darkmeson_mass_;
0141 
0142   // now do the angles, isotropic
0143   phi = 2.0 * M_PI * (fRndmPtr_->flat());
0144   theta = acos(2.0 * (fRndmPtr_->flat()) - 1.0);
0145 
0146   // compose the 4 std::vector
0147   en = sqrt(momentum * momentum + darkmeson_mass_ * darkmeson_mass_);
0148   Pythia8::Vec4 daughterFourMomentum =
0149       Pythia8::Vec4(momentum * cos(phi) * sin(theta), momentum * sin(phi) * sin(theta), momentum * sin(theta), en);
0150   return daughterFourMomentum;
0151 }
0152 
0153 // Auxiliary function, to be solved in order to impose energy conservation
0154 // To ballance energy, we solve for "scale" by demanding that this function vanishes, i.e. that shower energy and mediator energy are equal
0155 // By rescaling the momentum rather than the energy, avoid having to do these annoying rotations from the previous version
0156 const double SuepShower::reballanceFunction(double scale, const std::vector<Pythia8::Vec4>& shower) {
0157   double showerEnergy = 0.0;
0158   for (const auto& daughter : shower) {
0159     showerEnergy += sqrt(scale * scale * daughter.pAbs2() + darkmeson_mass_ * darkmeson_mass_);
0160   }
0161   return showerEnergy - mediator_energy_;
0162 }