No OneTemporary
Actions

Size

25 KB

Subscribers

None

View Options

	diff --git a/PDF/LHAPDF6.cc b/PDF/LHAPDF6.cc
	--- a/PDF/LHAPDF6.cc
	+++ b/PDF/LHAPDF6.cc
	@@ -1,512 +1,546 @@
	// -- C++ --
	//
	// LHAPDF6.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// // Copyright (C) 2014-2019 Leif Lonnblad, David Grellscheid
	//
	// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the LHAPDF class.
	//

	#include "LHAPDF6.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Command.h"
	#include "ThePEG/Interface/Deleted.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "LHAPDF/LHAPDF.h"

	using std::vector;
	using std::string;
	using std::pair;

	using ThePEG::GeV2;
	using ThePEG::cPDVector;

	ThePEG::LHAPDF::LHAPDF()
	: thePDF(), thePDFName("THEPEG_NO_PDFSET_CHOSEN"),
	theMember(0), theMaxFlav(5),
	xMin(0.0), xMax(1.0), Q2Min(ZERO), Q2Max(Constants::MaxEnergy2) {}

	ThePEG::LHAPDF::LHAPDF(const LHAPDF & x)
	: PDFBase(x),
	thePDF(), thePDFName(x.thePDFName),
	theMember(x.theMember),
	theMaxFlav(x.theMaxFlav),
	xMin(x.xMin), xMax(x.xMax), Q2Min(x.Q2Min), Q2Max(x.Q2Max) {}

	ThePEG::IBPtr ThePEG::LHAPDF::clone() const {
	return new_ptr(*this);
	}

	ThePEG::IBPtr ThePEG::LHAPDF::fullclone() const {
	return new_ptr(*this);
	}

	-void ThePEG::LHAPDF::initPDFptr() {
	+void ThePEG::LHAPDF::initPDFptr() const {
	::LHAPDF::setVerbosity(std::max(0,
	ThePEG::Debug::level - 1));
	if ( thePDF
	&& thePDF->set().name() == thePDFName
	&& thePDF->memberID() == theMember )
	return;
	delete thePDF;
	thePDF = ::LHAPDF::mkPDF(thePDFName, theMember);
	xMin = thePDF->xMin();
	xMax = thePDF->xMax();
	Q2Min = thePDF->q2Min() * GeV2;
	Q2Max = thePDF->q2Max() * GeV2;
	+ // get the PDG code for the particle for which this is the PDF
	+ thePID = thePDF->lhapdfID() < 10000 ? -211 : 2212;
	+ if(thePDF->info().has_key("Particle"))
	+ thePID = std::stoi(thePDF->info().get_entry("Particle"));
	}

	void ThePEG::LHAPDF::doinit() {
	PDFBase::doinit();
	initPDFptr();
	}

	void ThePEG::LHAPDF::dofinish() {
	PDFBase::dofinish();
	delete thePDF;
	thePDF = 0;
	}

	void ThePEG::LHAPDF::doinitrun() {
	PDFBase::doinitrun();
	initPDFptr();
	}

	void ThePEG::LHAPDF::setPDFName(string name) {
	if ( ::LHAPDF::endswith(name, ".LHgrid") ) {
	name = name.substr(0, name.size() - 7);
	}
	else if ( ::LHAPDF::endswith(name, ".LHpdf") ) {
	name = name.substr(0, name.size() - 6);
	}

	// fix the eternal typo
	if ( name == "cteq6ll" ) name = "cteq6l1";

	if ( ::LHAPDF::contains(::LHAPDF::availablePDFSets(), name) ) {
	thePDFName = name;
	theMember = 0;
	}
	else {
	Throw<ThePEG::LHAPDF::NotInstalled>()
	<< "'set " << fullName() << ":PDFName "
	<< name << "': PDF not installed. Try 'lhapdf install'.\n"
	<< Exception::setuperror;
	}
	}

	void ThePEG::LHAPDF::setPDFMember(int member) {
	try {
	::LHAPDF::PDFInfo * test =
	::LHAPDF::mkPDFInfo(thePDFName, member);
	if ( test )
	theMember = member;
	delete test;
	}
	catch (::LHAPDF::ReadError & e) {
	Throw<ThePEG::LHAPDF::NotInstalled>()
	<< e.what() << Exception::setuperror;
	}
	}

	string ThePEG::LHAPDF::doTest(string input) {
	double x = 0;
	Energy2 Q2 = ZERO;
	Energy2 P2 = ZERO;
	istringstream is(input);
	is >> x >> iunit(Q2, GeV2) >> iunit(P2, GeV2);
	initPDFptr();
	ostringstream os;
	for ( int i = 0; i < 13; ++i ) os << " " << thePDF->xfxQ2(i,x,Q2/GeV2);
	return os.str();
	}

	bool ThePEG::LHAPDF::canHandleParticle(tcPDPtr particle) const {
	+ initPDFptr();
	using namespace ParticleID;
	- return abs(particle->id()) == pplus \|\| abs(particle->id()) == n0;
	+ if(abs(thePID)==pplus)
	+ return abs(particle->id()) == pplus \|\| abs(particle->id()) == n0;
	+ else if(abs(thePID)==piplus)
	+ return abs(particle->id()) == piplus;
	+ else
	+ return false;
	}

	cPDVector ThePEG::LHAPDF::partons(tcPDPtr particle) const {
	// We assume that all standard partons can be extracted.
	const ::LHAPDF::PDFSet & pdfset = ::LHAPDF::getPDFSet(thePDFName);
	const vector<int> & flavs =
	pdfset.get_entry_as< vector<int> >("Flavors");

	// check for leptons and antileptons
	bool hasLepton [3] = {false,false,false};
	bool hasAntiLepton[3] = {false,false,false};
	vector<int> flav2;
	flav2.reserve(flavs.size()+3);
	for (int val : flavs) {
	if(abs(val) == 11 \|\| abs(val) ==13 \|\| abs(val)==15) {
	int ilep = (abs(val)-11)/2;
	if(val>0) hasLepton[ilep]=true;
	else hasAntiLepton[ilep]=true;
	}
	flav2.push_back(val);
	}
	for(unsigned int ix=0;ix<3;++ix) {
	if(hasLepton[ix] && !hasAntiLepton[ix]) {
	int pid = 11+2*ix;
	cerr << "Charged lepton with pid=" << pid << " found in PDF but not the antilepton"
	<< " assumming antilepton PDF is equal to the lepton one\n";
	flav2.push_back(-pid);
	}
	}

	cPDVector ret;
	ret.reserve( flav2.size() );
	if ( canHandleParticle(particle) ) {
	for ( size_t i=0; i < flav2.size(); ++i ) {
	ret.push_back(getParticleData(flav2[i]));
	}
	}
	assert( ret.size() == flav2.size() );
	return ret;
	}

	double ThePEG::LHAPDF::xfx(tcPDPtr particle, tcPDPtr parton,
	Energy2 partonScale,
	double x, double, Energy2) const {
	// Here we should return the actual density.
	using namespace ThePEG::ParticleID;

	double Q2 = partonScale/GeV2;

	if ( ! thePDF->inRangeXQ2(x, Q2) ) {
	switch ( rangeException ) {
	case rangeThrow: Throw<Exception>()
	<< "Momentum fraction (x=" << x << ") or scale (Q2=" << Q2
	<< " GeV^2) was outside of limits in PDF " << name() << "."
	<< Exception::eventerror;
	break;
	case rangeZero:
	return 0.0;
	case rangeFreeze:
	x = min(max(x, xMin), xMax);
	Q2 = min(max(Q2, Q2Min/GeV2), Q2Max/GeV2);
	}
	}

	int pid = parton->id();
	int abspid = abs(pid);

	switch ( pid ) {
	case t:
	case tbar:
	case b:
	case bbar:
	case c:
	case cbar:
	return maxFlav() < abspid ? 0.0 : thePDF->xfxQ2(pid,x,Q2);
	case s:
	case sbar:
	return thePDF->xfxQ2(pid,x,Q2);
	case u:
	switch ( particle->id() ) {
	case n0: return thePDF->xfxQ2(d ,x,Q2);
	case pbarminus: return thePDF->xfxQ2(ubar,x,Q2);
	case nbar0: return thePDF->xfxQ2(dbar,x,Q2);
	+ case piplus:
	+ case piminus: return thePDF->xfxQ2(u*thePID/particle->id(),x,Q2);
	case pplus:
	default: return thePDF->xfxQ2(u ,x,Q2);
	}
	case ubar:
	switch ( particle->id() ) {
	case n0: return thePDF->xfxQ2(dbar,x,Q2);
	case pbarminus: return thePDF->xfxQ2(u ,x,Q2);
	case nbar0: return thePDF->xfxQ2(d ,x,Q2);
	+ case piplus:
	+ case piminus: return thePDF->xfxQ2(ubar*thePID/particle->id(),x,Q2);
	case pplus:
	default: return thePDF->xfxQ2(ubar,x,Q2);
	}
	case d:
	switch ( particle->id() ) {
	case n0: return thePDF->xfxQ2(u ,x,Q2);
	case pbarminus: return thePDF->xfxQ2(dbar,x,Q2);
	case nbar0: return thePDF->xfxQ2(ubar,x,Q2);
	+ case piplus:
	+ case piminus: return thePDF->xfxQ2(d*thePID/particle->id(),x,Q2);
	case pplus:
	default: return thePDF->xfxQ2(d ,x,Q2);
	}
	case dbar:
	switch ( particle->id() ) {
	case n0: return thePDF->xfxQ2(ubar,x,Q2);
	case pbarminus: return thePDF->xfxQ2(d ,x,Q2);
	case nbar0: return thePDF->xfxQ2(u ,x,Q2);
	+ case piplus:
	+ case piminus: return thePDF->xfxQ2(dbar*thePID/particle->id(),x,Q2);
	case pplus:
	default: return thePDF->xfxQ2(dbar,x,Q2);
	}
	case g:
	return thePDF->xfxQ2(g,x,Q2);
	case ParticleID::gamma:
	return thePDF->xfxQ2(ParticleID::gamma,x,Q2);
	case ParticleID::eminus: case ParticleID::muminus: case ParticleID::tauminus:
	case ParticleID::eplus : case ParticleID::muplus : case ParticleID::tauplus :
	if(pid>0 \|\| thePDF->hasFlavor(pid))
	return thePDF->xfxQ2(pid,x,Q2);
	else
	return thePDF->xfxQ2(-pid,x,Q2);
	}
	return 0.0;
	}

	double ThePEG::LHAPDF::xfvl(tcPDPtr particle, tcPDPtr parton,
	Energy2 partonScale,
	double l, Energy2 particleScale) const {
	using Math::exp1m;
	return xfvx(particle, parton, partonScale,
	exp(-l), exp1m(-l), particleScale);
	}

	double ThePEG::LHAPDF::xfvx(tcPDPtr particle, tcPDPtr parton,
	Energy2 partonScale,
	double x, double, Energy2) const {
	// Here we should return the actual valence density. This will only
	// work properly for nucleons
	using namespace ThePEG::ParticleID;

	double Q2 = partonScale / GeV2;

	if ( ! thePDF->inRangeXQ2(x, Q2) ) {
	switch ( rangeException ) {
	case rangeThrow: Throw<Exception>()
	<< "Momentum fraction (x=" << x << ") or scale (Q2=" << Q2
	<< " GeV^2) was outside of limits in PDF " << name() << "."
	<< Exception::eventerror;
	break;
	case rangeZero:
	return 0.0;
	case rangeFreeze:
	x = min(max(x, xMin), xMax);
	Q2 = min(max(Q2, Q2Min/GeV2), Q2Max/GeV2);
	}
	}

	switch ( parton->id() ) {
	case t:
	case tbar:
	case b:
	case bbar:
	case c:
	case cbar:
	case s:
	case sbar:
	case ParticleID::gamma:
	return 0.0;
	case u:
	switch ( particle->id() ) {
	case n0: return thePDF->xfxQ2(d,x,Q2) - thePDF->xfxQ2(dbar,x,Q2);
	case pbarminus: return 0.0;
	case nbar0: return 0.0;
	case pplus: return thePDF->xfxQ2(u,x,Q2) - thePDF->xfxQ2(ubar,x,Q2);
	+ case piplus: return thePDF->xfxQ2(thePID/particle->id()u,x,Q2)-thePDF->xfxQ2(thePID/particle->id()ubar,x,Q2);
	+ case piminus: return 0.;
	default: return 0.0;
	}
	case ubar:
	switch ( particle->id() ) {
	case n0: return 0.0;
	case pbarminus: return thePDF->xfxQ2(u,x,Q2) - thePDF->xfxQ2(ubar,x,Q2);
	case nbar0: return thePDF->xfxQ2(d,x,Q2) - thePDF->xfxQ2(dbar,x,Q2);
	+ case piplus: return 0.;
	+ case piminus: return thePDF->xfxQ2(thePID/particle->id()ubar,x,Q2)-thePDF->xfxQ2(thePID/particle->id()u,x,Q2);
	case pplus:
	default: return 0.0;
	}
	case d:
	switch ( particle->id() ) {
	case n0: return thePDF->xfxQ2(u,x,Q2) - thePDF->xfxQ2(ubar,x,Q2);
	case pbarminus: return 0.0;
	case nbar0: return 0.0;
	case pplus: return thePDF->xfxQ2(d,x,Q2) - thePDF->xfxQ2(dbar,x,Q2);
	+ case piplus: return 0.;
	+ case piminus: return thePDF->xfxQ2(thePID/particle->id()d,x,Q2)-thePDF->xfxQ2(thePID/particle->id()dbar,x,Q2);
	default: return 0.0;
	}
	case dbar:
	switch ( particle->id() ) {
	case n0: return 0.0;
	case pbarminus: return thePDF->xfxQ2(d,x,Q2) - thePDF->xfxQ2(dbar,x,Q2);
	case nbar0: return thePDF->xfxQ2(u,x,Q2) - thePDF->xfxQ2(ubar,x,Q2);
	+ case piplus: return thePDF->xfxQ2(thePID/particle->id()dbar,x,Q2)-thePDF->xfxQ2(thePID/particle->id()d,x,Q2);
	+ case piminus: return 0.;
	case pplus:
	default: return 0.0;
	}
	case ParticleID::g:
	return 0.0;
	}
	return 0.0;
	}

	double ThePEG::LHAPDF::xfsx(tcPDPtr particle, tcPDPtr parton,
	Energy2 partonScale, double x,
	double, Energy2) const {
	// Here we should return the actual density.
	using namespace ThePEG::ParticleID;

	double Q2 = partonScale / GeV2;

	if ( ! thePDF->inRangeXQ2(x, Q2) ) {
	switch ( rangeException ) {
	case rangeThrow: Throw<Exception>()
	<< "Momentum fraction (x=" << x << ") or scale (Q2=" << Q2
	<< " GeV^2) was outside of limits in PDF " << name() << "."
	<< Exception::eventerror;
	break;
	case rangeZero:
	return 0.0;
	case rangeFreeze:
	x = min(max(x, xMin), xMax);
	Q2 = min(max(Q2, Q2Min/GeV2), Q2Max/GeV2);
	}
	}
	int pid = parton->id();
	int abspid = abs(pid);

	switch ( pid ) {
	case t:
	case tbar:
	case b:
	case bbar:
	case c:
	case cbar:
	return maxFlav() < abspid ? 0.0 : thePDF->xfxQ2(pid,x,Q2);
	case s:
	case sbar:
	return thePDF->xfxQ2(pid,x,Q2);
	case u:
	switch ( particle->id() ) {
	case n0: return thePDF->xfxQ2(dbar,x,Q2);
	case pbarminus: return thePDF->xfxQ2(ubar,x,Q2);
	case nbar0: return thePDF->xfxQ2(dbar,x,Q2);
	case pplus: return thePDF->xfxQ2(ubar,x,Q2);
	+ case piplus: return thePDF->xfxQ2(thePID/particle->id()*ubar,x,Q2);
	+ case piminus: return thePDF->xfxQ2(thePID/particle->id()*u ,x,Q2);
	default: return thePDF->xfxQ2(u ,x,Q2);
	}
	case ubar:
	switch ( particle->id() ) {
	case n0: return thePDF->xfxQ2(dbar,x,Q2);
	case pbarminus: return thePDF->xfxQ2(ubar,x,Q2);
	case nbar0: return thePDF->xfxQ2(dbar,x,Q2);
	+ case piplus: return thePDF->xfxQ2(thePID/particle->id()*ubar,x,Q2);
	+ case piminus: return thePDF->xfxQ2(thePID/particle->id()*u ,x,Q2);
	case pplus:
	default: return thePDF->xfxQ2(ubar,x,Q2);
	}
	case d:
	switch ( particle->id() ) {
	case n0: return thePDF->xfxQ2(ubar,x,Q2);
	case pbarminus: return thePDF->xfxQ2(dbar,x,Q2);
	case nbar0: return thePDF->xfxQ2(ubar,x,Q2);
	case pplus: return thePDF->xfxQ2(dbar,x,Q2);
	+ case piplus: return thePDF->xfxQ2(thePID/particle->id()*d ,x,Q2);
	+ case piminus: return thePDF->xfxQ2(thePID/particle->id()*dbar,x,Q2);
	default: return thePDF->xfxQ2(d ,x,Q2);
	}
	case dbar:
	switch ( particle->id() ) {
	case n0: return thePDF->xfxQ2(ubar,x,Q2);
	case pbarminus: return thePDF->xfxQ2(dbar,x,Q2);
	case nbar0: return thePDF->xfxQ2(ubar,x,Q2);
	+ case piplus: return thePDF->xfxQ2(thePID/particle->id()*d ,x,Q2);
	+ case piminus: return thePDF->xfxQ2(thePID/particle->id()*dbar,x,Q2);
	case pplus:
	default: return thePDF->xfxQ2(dbar,x,Q2);
	}
	case g:
	return thePDF->xfxQ2(g,x,Q2);
	case ParticleID::gamma:
	return thePDF->xfxQ2(ParticleID::gamma,x,Q2);
	case ParticleID::eminus: case ParticleID::muminus: case ParticleID::tauminus:
	case ParticleID::eplus : case ParticleID::muplus : case ParticleID::tauplus :
	if(pid>0 \|\| thePDF->hasFlavor(pid))
	return thePDF->xfxQ2(pid,x,Q2);
	else
	return thePDF->xfxQ2(-pid,x,Q2);
	}
	return 0.0;
	}

	void ThePEG::LHAPDF::persistentOutput(PersistentOStream & os) const {
	os << thePDFName << theMember << theMaxFlav
	<< xMin << xMax << ounit(Q2Min, GeV2) << ounit(Q2Max, GeV2);
	}

	void ThePEG::LHAPDF::persistentInput(PersistentIStream & is, int) {
	is >> thePDFName >> theMember >> theMaxFlav
	>> xMin >> xMax >> iunit(Q2Min, GeV2) >> iunit(Q2Max, GeV2);
	initPDFptr();
	}

	ThePEG::DescribeClass<ThePEG::LHAPDF, ThePEG::PDFBase>
	describeLHAPDF("ThePEG::LHAPDF", "ThePEGLHAPDF.so");

	void ThePEG::LHAPDF::Init() {

	static ClassDocumentation<LHAPDF> documentation
	("The LHAPDF class inherits from PDFBase and implements an interface "
	"to the LHAPDF library of parton density function parameterizations. "
	"This class is available even if LHAPDF was not properly installed "
	"when ThePEG was installed, but will then produce an error in the "
	"initialization. Note that the valence densities from the xfvx() and "
	"xfvl() function will only work properly for nucleons. All other "
	"particles will have zero valence densities.");

	static Deleted<LHAPDF> interfacePType
	("PType",
	"The LHAPDFv6 interface currently does not support pi.");

	static Parameter<LHAPDF,string> interfacePDFName
	("PDFName",
	"The name of the PDF set to be used. Should correspond to "
	"the LHAPDF v6 name.",
	&ThePEG::LHAPDF::thePDFName, "THEPEG_NO_PDFSET_CHOSEN", true, false,
	&ThePEG::LHAPDF::setPDFName);

	static Deleted<LHAPDF> interfacePDFNumber
	("PDFNumber",
	"Not implemented in the LHAPDFv6 interface. "
	"Use :PDFName and :Member instead.");

	static Parameter<LHAPDF,int> interfaceMember
	("Member",
	"The chosen member of the selected PDF set.",
	&ThePEG::LHAPDF::theMember, 0, 0, 0,
	true, false, Interface::lowerlim,
	&ThePEG::LHAPDF::setPDFMember);

	static Deleted<LHAPDF> interfacePDFLIBNumbers
	("PDFLIBNumbers",
	"Not implemented in the LHAPDFv6 interface. "
	"Use :PDFName and :Member instead.");

	static Deleted<LHAPDF> interfaceEnablePartonicGamma
	("EnablePartonicGamma",
	"Not required in LHAPDFv6.");

	static Deleted<LHAPDF> interfacePhotonOption
	("PhotonOption",
	"Not required in LHAPDFv6.");

	static Deleted<LHAPDF> interfaceMaxNSet
	("MaxNSet",
	"Not required in LHAPDFv6.");

	static Command<LHAPDF> interfaceTest
	("Test",
	"Write out the values of the chosen PDF set using the x, Q2 and P2 "
	"parameters supplied.",
	&ThePEG::LHAPDF::doTest, true);

	static Deleted<LHAPDF> interfaceVerboseLevel
	("VerboseLevel",
	"LHAPDFv6 uses general debug level instead.");

	static Parameter<LHAPDF,int> interfaceMaxFlav
	("MaxFlav",
	"The maximum number of flavours for which non-zero densities are "
	"reported. The actual number of flavours may be less depending on "
	"the chosen PDF set.",
	&ThePEG::LHAPDF::theMaxFlav, 5, 3, 0,
	true, false, Interface::lowerlim);

	interfacePDFName.rank(9);
	interfaceMember.rank(8);

	}

	diff --git a/PDF/LHAPDF6.h b/PDF/LHAPDF6.h
	--- a/PDF/LHAPDF6.h
	+++ b/PDF/LHAPDF6.h
	@@ -1,299 +1,301 @@
	// -- C++ --
	//
	// LHAPDF6.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 2014-2019 Leif Lonnblad, David Grellscheid
	//
	// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef THEPEG_LHAPDF6_H
	#define THEPEG_LHAPDF6_H
	//
	// This is the declaration of the LHAPDF class.
	//

	#include "ThePEG/PDF/PDFBase.h"

	namespace LHAPDF {
	- class PDF;
	+ class PDF;
	}

	namespace ThePEG {

	/**
	* The LHAPDF class inherits from PDFBase and implements an interface
	* to the LHAPDF library of parton density function
	* parameterizations. This class is available even if LHAPDF was not
	* properly installed when ThePEG was installed, but will then produce
	* an error in the initialization.
	*
	* Note that the valence densities from the xfvx() and xfvl() function
	* will only work properly for nucleons. All other particles will have
	* zero valence densities.
	*
	* @see \ref LHAPDFInterfaces "The interfaces"
	* defined for LHAPDF.
	*/
	class LHAPDF: public PDFBase {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	LHAPDF();

	/**
	* The copy constructor.
	*/
	LHAPDF(const LHAPDF &);
	//@}

	public:

	/** @name Virtual functions to be overridden by sub-classes. */
	//@{
	/**
	* Return true if this PDF can handle the extraction of partons from
	* the given \a particle.
	*/
	virtual bool canHandleParticle(tcPDPtr particle) const;

	/**
	* Return the partons which this PDF may extract from the given
	* \a particle.
	*/
	virtual cPDVector partons(tcPDPtr particle) const;

	/**
	* The density. Return the pdf for the given \a parton inside the
	* given \a particle for the virtuality \a partonScale and momentum
	* fraction \a x (with x = 1-\a eps). The \a particle is assumed to
	* have a virtuality \a particleScale.
	*/
	virtual double xfx(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale,
	double x, double eps = 0.0,
	Energy2 particleScale = ZERO) const;

	/**
	* The valence density. Return the pdf for the given cvalence \a
	* parton inside the given \a particle for the virtuality \a
	* partonScale and logarithmic momentum fraction \a l. The \a
	* particle is assumed to have a virtuality \a particleScale. This
	* will only work properly for nucleons. All other particles will
	* have zero valence densities
	*/
	virtual double xfvl(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale,
	double l, Energy2 particleScale = ZERO) const;

	/**
	* The valence density. Return the pdf for the given cvalence \a
	* parton inside the given \a particle for the virtuality \a
	* partonScale and momentum fraction \a x (with x = 1-\a eps). The
	* \a particle is assumed to have a virtuality \a
	* particleScale. This will only work properly for nucleons. All
	* other particles will have zero valence densities
	*/
	virtual double xfvx(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale,
	double x, double eps = 0.0,
	Energy2 particleScale = ZERO) const;

	/**
	* The sea density. Return the pdf for the given cvalence \a
	* parton inside the given \a particle for the virtuality \a
	* partonScale and momentum fraction \a x. The \a particle is
	* assumed to have a virtuality \a particleScale. If not overidden
	* by a sub class this implementation will assume that the
	* difference between a quark and anti-quark distribution is due do
	* valense quarks.
	*/
	virtual double xfsx(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale,
	double x, double eps = 0.0,
	Energy2 particleScale = ZERO) const;
	//@}


	/** @name Simple access function. */
	//@{
	/**
	* The name if the PDF set to be used. The full name including the
	* <code>.LHpdf</code> or <code>.LHgrid</code> suffix.
	*/
	const string & PDFName() const { return thePDFName; }

	/**
	* The chosen member of the selected PDF set.
	*/
	int member() const { return theMember; }

	/**
	* The maximum number of flavours for which non-zero densities are
	* reported. The actual number of flavours may be less depending on
	* the chosen PDF set.
	*/
	int maxFlav() const { return theMaxFlav; }
	//@}

	protected:

	/** @name Internal helper functions. */
	//@{
	/**
	* Initialize the LHAPDF library for the chosen PDF set if it has
	* not been done before.
	*/
	- void initPDFptr();
	+ void initPDFptr() const;

	/**
	* Used by the interface to select a set according to a file name.
	*/
	void setPDFName(string name);

	/**
	* Used by the interface to select a member in the current set.
	*/
	void setPDFMember(int n);

	/**
	* Interface for simple tests.
	*/
	string doTest(string input);
	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	public:

	/** @cond EXCEPTIONCLASSES */
	/** Exception class used if the LHAPDF library was not installed. */
	class NotInstalled: public InterfaceException {};

	/** @endcond */

	-// If needed, insert declarations of virtual function defined in the
	-// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
	-
	protected:

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();

	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();

	/**
	* Finalize this object. Called in the run phase just after a
	* run has ended. Used eg. to write out statistics.
	*/
	virtual void dofinish();
	//@}

	private:

	/**
	* LHAPDF member object
	*/
	- ::LHAPDF::PDF * thePDF;
	+ mutable ::LHAPDF::PDF * thePDF;

	/**
	* The name of the selected PDF set.
	*/
	string thePDFName;

	/**
	* The chosen member of the selected PDF set.
	*/
	int theMember;

	/**
	* The maximum number of flavours for which non-zero densities are
	* reported. The actual number of flavours may be less depending on
	* the chosen PDF set.
	*/
	int theMaxFlav;

	/**
	+ * The PDG code for the particle described by the PDF
	+ */
	+ mutable long thePID;
	+
	+ /**
	* The minimum \f$x\f$-value for the current PDF set.
	*/
	- double xMin;
	+ mutable double xMin;

	/**
	* The maximum \f$x\f$-value for the current PDF set.
	*/
	- double xMax;
	+ mutable double xMax;

	/**
	* The minimum \f$Q^2\f$-value for the current PDF set.
	*/
	- Energy2 Q2Min;
	+ mutable Energy2 Q2Min;

	/**
	* The maximum \f$Q^2\f$-value for the current PDF set.
	*/
	- Energy2 Q2Max;
	+ mutable Energy2 Q2Max;

	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	LHAPDF & operator=(const LHAPDF &) = delete;

	};

	}

	#endif /* THEPEG_LHAPDF6_H */

File Metadata

Mime Type: text/x-diff
Expires: Fri, Apr 4, 9:38 PM (16 h, 2 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4737426
Default Alt Text: (25 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions