SuepShower Class Reference

#include <SuepShower.h>

Public Member Functions
std::vector< Pythia8::Vec4 >	generateShower (double energy)
	SuepShower (double mass, double temperature, Pythia8::Rndm *rndmPtr)
	~SuepShower ()

Private Member Functions
const double	fMaxwellBoltzmann (double p)
const double	fMaxwellBoltzmannPrime (double p)
const Pythia8::Vec4	generateFourVector ()
const double	logTestFunction (double p)
const double	reballanceFunction (double scale, const std::vector< Pythia8::Vec4 > &shower)

Private Attributes
double	darkmeson_mass_
Pythia8::Rndm *	fRndmPtr_
double	lambda_minus_
double	lambda_plus_
double	mass_over_T_
double	mediator_energy_
double	p_m_
double	p_minus_
double	p_plus_
double	q_m_
double	q_minus_
double	q_plus_
boost::math::tools::eps_tolerance< double >	tolerance_

Detailed Description

Definition at line 12 of file SuepShower.h.

Constructor & Destructor Documentation

SuepShower()

SuepShower::SuepShower	(	double	mass,
		double	temperature,
		Pythia8::Rndm *	rndmPtr
	)

Definition at line 18 of file SuepShower.cc.

References darkmeson_mass_, first, fMaxwellBoltzmann(), fMaxwellBoltzmannPrime(), fRndmPtr_, lambda_minus_, lambda_plus_, logTestFunction(), EgHLTOffHistBins_cfi::mass, mass_over_T_, p_m_, p_minus_, p_plus_, q_m_, q_minus_, q_plus_, mathSSE::sqrt(), jvcParameters_cfi::temperature, and tolerance_.

                                                                            {
  edm::LogInfo("Pythia8Interface") << "Creating a SuepShower module";
  // Parameter setup
  darkmeson_mass_ = mass;
  fRndmPtr_ = rndmPtr;
  // No need to save temperature, as everything depends on it through m/T
  mass_over_T_ = darkmeson_mass_ / temperature;

  // 128 bit numerical precision, i.e. take machine precision
  tolerance_ = boost::math::tools::eps_tolerance<double>(128);

  // Median momentum
  p_m_ = sqrt(2 / (mass_over_T_ * mass_over_T_) * (1 + sqrt(1 + mass_over_T_ * mass_over_T_)));
  double pmax = sqrt(2 + 2 * sqrt(1 + mass_over_T_ * mass_over_T_)) / mass_over_T_;
  // Find the two cases where f(p)/f(pmax) = 1/e, given MB distribution, one appears at each side of pmax always
  p_plus_ = (boost::math::tools::bisect(
                 boost::bind(&SuepShower::logTestFunction, this, boost::placeholders::_1), pmax, 50.0, tolerance_))
                .first;  // first root
  p_minus_ = (boost::math::tools::bisect(
                  boost::bind(&SuepShower::logTestFunction, this, boost::placeholders::_1), 0.0, pmax, tolerance_))
                 .first;  // second root
  // Define the auxiliar quantities for the random generation of momenta
  lambda_plus_ = -fMaxwellBoltzmann(p_plus_) / fMaxwellBoltzmannPrime(p_plus_);
  lambda_minus_ = fMaxwellBoltzmann(p_minus_) / fMaxwellBoltzmannPrime(p_minus_);
  q_plus_ = lambda_plus_ / (p_plus_ - p_minus_);
  q_minus_ = lambda_minus_ / (p_plus_ - p_minus_);
  q_m_ = 1 - (q_plus_ + q_minus_);
}

~SuepShower()

SuepShower::~SuepShower

(

)

Definition at line 47 of file SuepShower.cc.

{}

Member Function Documentation

fMaxwellBoltzmann()

const double SuepShower::fMaxwellBoltzmann ( double  p )

private

Definition at line 102 of file SuepShower.cc.

References JetChargeProducer_cfi::exp, mass_over_T_, AlCaHLTBitMon_ParallelJobs::p, and mathSSE::sqrt().

Referenced by generateFourVector(), logTestFunction(), and SuepShower().

                                                   {
  return p * p * exp(-mass_over_T_ * p * p / (1 + sqrt(1 + p * p)));
}

fMaxwellBoltzmannPrime()

const double SuepShower::fMaxwellBoltzmannPrime ( double  p )

private

Definition at line 107 of file SuepShower.cc.

References JetChargeProducer_cfi::exp, mass_over_T_, AlCaHLTBitMon_ParallelJobs::p, and mathSSE::sqrt().

Referenced by SuepShower().

                                                        {
  return exp(-mass_over_T_ * p * p / (1 + sqrt(1 + p * p))) * p * (2 - mass_over_T_ * p * p / sqrt(1 + p * p));
}

generateFourVector()

const Pythia8::Vec4 SuepShower::generateFourVector ( )

private

Definition at line 115 of file SuepShower.cc.

References funct::cos(), darkmeson_mass_, JetChargeProducer_cfi::exp, fMaxwellBoltzmann(), fRndmPtr_, mps_fire::i, lambda_minus_, lambda_plus_, dqm-mbProfile::log, M_PI, p_m_, p_minus_, p_plus_, phi, q_m_, q_minus_, q_plus_, funct::sin(), mathSSE::sqrt(), theta(), mitigatedMETSequence_cff::U, cms::cuda::V, X, and BeamSpotPI::Y.

Referenced by generateShower().

                                                 {
  double en, phi, theta, momentum;  //kinematic variables of the 4 std::vector

  // First do momentum, following naming of arxiv:1305.5226
  double U, V, X, Y, E;
  int i = 0;
  while (i < 100) {
    // Very rarely (<1e-5 events) takes more than one loop to converge, set to 100 for safety
    U = fRndmPtr_->flat();
    V = fRndmPtr_->flat();

    if (U < q_m_) {
      Y = U / q_m_;
      X = (1 - Y) * (p_minus_ + lambda_minus_) + Y * (p_plus_ - lambda_plus_);
      if (V < fMaxwellBoltzmann(X) / fMaxwellBoltzmann(p_m_) && X > 0) {
        break;
      }
    } else {
      if (U < q_m_ + q_plus_) {
        E = -log((U - q_m_) / q_plus_);
        X = p_plus_ - lambda_plus_ * (1 - E);
        if (V < exp(E) * fMaxwellBoltzmann(X) / fMaxwellBoltzmann(p_m_) && X > 0) {
          break;
        }
      } else {
        E = -log((U - (q_m_ + q_plus_)) / q_minus_);
        X = p_minus_ + lambda_minus_ * (1 - E);
        if (V < exp(E) * fMaxwellBoltzmann(X) / fMaxwellBoltzmann(p_m_) && X > 0) {
          break;
        }
      }
    }
  }
  // X is the dimensionless momentum, p/m
  momentum = X * darkmeson_mass_;

  // now do the angles, isotropic
  phi = 2.0 * M_PI * (fRndmPtr_->flat());
  theta = acos(2.0 * (fRndmPtr_->flat()) - 1.0);

  // compose the 4 std::vector
  en = sqrt(momentum * momentum + darkmeson_mass_ * darkmeson_mass_);
  Pythia8::Vec4 daughterFourMomentum =
      Pythia8::Vec4(momentum * cos(phi) * sin(theta), momentum * sin(phi) * sin(theta), momentum * sin(theta), en);
  return daughterFourMomentum;
}

generateShower()

std::vector< Pythia8::Vec4 > SuepShower::generateShower

(

double

energy

)

Definition at line 50 of file SuepShower.cc.

References pfMETCorrectionType0_cfi::correction, darkmeson_mass_, MillePedeFileConverter_cfg::e, HCALHighEnergyHPDFilter_cfi::energy, first, generateFourVector(), mediator_energy_, reballanceFunction(), L1EGammaClusterEmuProducer_cfi::scale, mathSSE::sqrt(), and tolerance_.

                                                               {
  mediator_energy_ = energy;
  std::vector<Pythia8::Vec4> shower;
  double shower_energy = 0.0;

  // Fill up shower record
  while (shower_energy < mediator_energy_ || shower.size() < 2) {
    shower.push_back(generateFourVector());
    shower_energy += (shower.back()).e();
  }

  // Reballance momenta to ensure conservation
  // Correction starts at 0
  Pythia8::Vec4 correction = Pythia8::Vec4(0., 0., 0., 0.);
  for (const auto& daughter : shower) {
    correction = correction + daughter;
  }
  correction = correction / shower.size();
  // We only want to correct momenta first
  correction.e(0);
  for (auto& daughter : shower) {
    daughter = daughter - correction;
  }

  // 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
  double scale;
  double minscale = 0.0;
  double maxscale = 2.0;
  while (SuepShower::reballanceFunction(minscale, shower) * SuepShower::reballanceFunction(maxscale, shower) > 0) {
    minscale = maxscale;
    maxscale *= 2;
  }

  scale =
      (boost::math::tools::bisect(boost::bind(&SuepShower::reballanceFunction, this, boost::placeholders::_1, shower),
                                  minscale,
                                  maxscale,
                                  tolerance_))
          .first;

  for (auto& daughter : shower) {
    daughter.px(daughter.px() * scale);
    daughter.py(daughter.py() * scale);
    daughter.pz(daughter.pz() * scale);
    // Force everything to be on-shell
    daughter.e(sqrt(daughter.pAbs2() + darkmeson_mass_ * darkmeson_mass_));
  }
  return shower;
}

logTestFunction()

const double SuepShower::logTestFunction ( double  p )

private

Definition at line 112 of file SuepShower.cc.

References fMaxwellBoltzmann(), dqm-mbProfile::log, AlCaHLTBitMon_ParallelJobs::p, and p_m_.

Referenced by SuepShower().

{ return log(fMaxwellBoltzmann(p) / fMaxwellBoltzmann(p_m_)) + 1.0; }

reballanceFunction()

const double SuepShower::reballanceFunction ( double  scale, const std::vector< Pythia8::Vec4 > &  shower  )

private

Definition at line 165 of file SuepShower.cc.

References darkmeson_mass_, mediator_energy_, L1EGammaClusterEmuProducer_cfi::scale, and mathSSE::sqrt().

Referenced by generateShower().

                                                                                              {
  double showerEnergy = 0.0;
  for (const auto& daughter : shower) {
    showerEnergy += sqrt(scale * scale * daughter.pAbs2() + darkmeson_mass_ * darkmeson_mass_);
  }
  return showerEnergy - mediator_energy_;
}

Member Data Documentation

darkmeson_mass_

double SuepShower::darkmeson_mass_

private

Definition at line 26 of file SuepShower.h.

Referenced by generateFourVector(), generateShower(), reballanceFunction(), and SuepShower().

fRndmPtr_

Pythia8::Rndm* SuepShower::fRndmPtr_

private

Definition at line 46 of file SuepShower.h.

Referenced by generateFourVector(), and SuepShower().

lambda_minus_

double SuepShower::lambda_minus_

private

Definition at line 41 of file SuepShower.h.

Referenced by generateFourVector(), and SuepShower().

lambda_plus_

double SuepShower::lambda_plus_

private

Definition at line 41 of file SuepShower.h.

Referenced by generateFourVector(), and SuepShower().

mass_over_T_

double SuepShower::mass_over_T_

private

Definition at line 27 of file SuepShower.h.

Referenced by fMaxwellBoltzmann(), fMaxwellBoltzmannPrime(), and SuepShower().

mediator_energy_

double SuepShower::mediator_energy_

private

Definition at line 30 of file SuepShower.h.

Referenced by generateShower(), and reballanceFunction().

p_m_

double SuepShower::p_m_

private

Definition at line 37 of file SuepShower.h.

Referenced by generateFourVector(), logTestFunction(), and SuepShower().

p_minus_

double SuepShower::p_minus_

private

Definition at line 39 of file SuepShower.h.

Referenced by generateFourVector(), and SuepShower().

p_plus_

double SuepShower::p_plus_

private

Definition at line 39 of file SuepShower.h.

Referenced by generateFourVector(), and SuepShower().

q_m_

double SuepShower::q_m_

private

Definition at line 43 of file SuepShower.h.

Referenced by generateFourVector(), and SuepShower().

q_minus_

double SuepShower::q_minus_

private

Definition at line 43 of file SuepShower.h.

Referenced by generateFourVector(), and SuepShower().

q_plus_

double SuepShower::q_plus_

private

Definition at line 43 of file SuepShower.h.

Referenced by generateFourVector(), and SuepShower().

tolerance_

boost::math::tools::eps_tolerance<double> SuepShower::tolerance_

private

Definition at line 33 of file SuepShower.h.

Referenced by generateShower(), and SuepShower().