#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	(	double	mass,
		double	temperature,
		Pythia8::Rndm *	rndmPtr
	)

Definition at line 18 of file SuepShower.cc.

References darkmeson_mass_, plotBeamSpotDB::first, fMaxwellBoltzmann(), fMaxwellBoltzmannPrime(), fRndmPtr_, lambda_minus_, lambda_plus_, logTestFunction(), ResonanceBuilder::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 ( )

Definition at line 47 of file SuepShower.cc.

47 {}

Member Function Documentation

const double SuepShower::fMaxwellBoltzmann ( double p )

private

Definition at line 93 of file SuepShower.cc.

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

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

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

const double SuepShower::fMaxwellBoltzmannPrime ( double p )

private

Definition at line 98 of file SuepShower.cc.

References JetChargeProducer_cfi::exp, mass_over_T_, 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));
 }

const Pythia8::Vec4 SuepShower::generateFourVector ( )

private

Definition at line 106 of file SuepShower.cc.

References funct::cos(), darkmeson_mass_, JetChargeProducer_cfi::exp, fMaxwellBoltzmann(), fRndmPtr_, mps_fire::i, lambda_minus_, lambda_plus_, cmsBatch::log, M_PI, p_m_, p_minus_, p_plus_, phi, q_m_, q_minus_, q_plus_, funct::sin(), mathSSE::sqrt(), theta(), mitigatedMETSequence_cff::U, X, and DOFs::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;
 }

std::vector< Pythia8::Vec4 > SuepShower::generateShower ( double energy )

Definition at line 50 of file SuepShower.cc.

References darkmeson_mass_, MillePedeFileConverter_cfg::e, plotBeamSpotDB::first, generateFourVector(), mediator_energy_, reballanceFunction(), Scenarios_cff::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.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;
   scale =
       (boost::math::tools::bisect(
            boost::bind(&SuepShower::reballanceFunction, this, boost::placeholders::_1, shower), 0.0, 2.0, 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;
 }

const double SuepShower::logTestFunction ( double p )

private

Definition at line 103 of file SuepShower.cc.

References fMaxwellBoltzmann(), cmsBatch::log, and p_m_.

Referenced by SuepShower().

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

AlCaHLTBitMon_ParallelJobs.p

p

Definition: AlCaHLTBitMon_ParallelJobs.py:153

SuepShower::p_m_

double p_m_

Definition: SuepShower.h:37

cmsBatch.log

log

Definition: cmsBatch.py:343

SuepShower::fMaxwellBoltzmann

const double fMaxwellBoltzmann(double p)

Definition: SuepShower.cc:93

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

private

Definition at line 156 of file SuepShower.cc.

References darkmeson_mass_, mediator_energy_, 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_;
 }