3
// This file is part of LHAPDF
4
// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
6
#include "LHAPDF/PDF.h"
7
#include "LHAPDF/PDFSet.h"
8
#include "LHAPDF/PDFIndex.h"
9
#include "LHAPDF/Factories.h"
10
#include "LHAPDF/Utils.h"
11
#include "LHAPDF/Paths.h"
12
#include "LHAPDF/Version.h"
13
#include "LHAPDF/LHAGlue.h"
18
// We have to create and initialise some common blocks here for backwards compatibility
25
double xmin, xmax, q2min, q2max;
30
double qcdlha4, qcdlha5;
37
namespace lhapdf_amc { //< Unnamed namespace to restrict visibility to this file
39
/// @brief PDF object storage here is a smart pointer to ensure deletion of created PDFs
41
/// NB. std::auto_ptr cannot be stored in STL containers, hence we use
42
/// boost::shared_ptr. std::unique_ptr is the nature replacement when C++11
43
/// is globally available.
44
typedef boost::shared_ptr<LHAPDF::PDF> PDFPtr;
46
/// @brief A struct for handling the active PDFs for the Fortran interface.
48
/// We operate in a string-based way, since maybe there will be sets with names, but no
49
/// index entry in pdfsets.index.
51
/// @todo Can we avoid the strings and just work via the LHAPDF ID and factory construction?
53
/// Smart pointers are used in the native map used for PDF member storage so
54
/// that they auto-delete if the PDFSetHandler that holds them goes out of
55
/// scope (i.e. is overwritten).
56
struct PDFSetHandler {
58
/// Default constructor
59
PDFSetHandler() : currentmem(0)
60
{ } //< It'll be stored in a map so we need one of these...
62
/// Constructor from a PDF set name
63
PDFSetHandler(const string& name)
69
/// Constructor from a PDF set's LHAPDF ID code
70
PDFSetHandler(int lhaid) {
71
pair<string,int> set_mem = LHAPDF::lookupPDF(lhaid);
72
// First check that the lookup was successful, i.e. it was a valid ID for the LHAPDF6 set collection
73
if (set_mem.first.empty() || set_mem.second < 0)
74
throw LHAPDF::UserError("Could not find a valid PDF with LHAPDF ID = " + LHAPDF::to_str(lhaid));
75
// Try to load this PDF (checking that the member number is in the set's range is done in mkPDF, called by loadMember)
76
setname = set_mem.first;
77
loadMember(set_mem.second);
80
/// @brief Load a new PDF member
82
/// If it's already loaded, the existing object will not be reloaded.
83
void loadMember(int mem) {
85
throw LHAPDF::UserError("Tried to load a negative PDF member ID: " + LHAPDF::to_str(mem) + " in set " + setname);
86
if (members.find(mem) == members.end())
87
members[mem] = PDFPtr(LHAPDF::mkPDF(setname, mem));
91
/// Actively delete a PDF member to save memory
92
void unloadMember(int mem) {
94
const int nextmem = (!members.empty()) ? members.begin()->first : 0;
98
/// @brief Get a PDF member
100
/// Non-const because it can secretly load the member. Not that constness
101
/// matters in a Fortran interface utility function!
102
const PDFPtr member(int mem) {
104
return members.find(mem)->second;
107
/// Get the currently active PDF member
109
/// Non-const because it can secretly load the member. Not that constness
110
/// matters in a Fortran interface utility function!
111
const PDFPtr activemember() {
112
return member(currentmem);
115
/// The currently active member in this set
121
/// Map of pointers to selected member PDFs
123
// /// It's mutable so that a "const" member-getting operation can implicitly
124
// /// load a new PDF object. Good idea / bad idea? Disabled for now.
125
// mutable map<int, PDFPtr> members;
126
map<int, PDFPtr> members;
130
/// Collection of active sets
131
static map<int, PDFSetHandler> ACTIVESETS;
133
/// The currently active set
140
string lhaglue_get_current_pdf(int nset) {
141
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
143
lhapdf_amc::CURRENTSET = nset;
144
return lhapdf_amc::ACTIVESETS[nset].activemember()->set().name() + " (" +
145
LHAPDF::to_str(lhapdf_amc::ACTIVESETS[nset].activemember()->lhapdfID()) + ")";
152
// NEW FORTRAN INTERFACE FUNCTIONS
154
/// List of available sets
155
void lhapdf_getversion_(char* s, size_t len) {
156
strncpy(s, LHAPDF_VERSION, len);
159
/// List of available PDF sets, returned as a space-separated string
160
void lhapdf_getpdfsetlist_(char* s, size_t len) {
162
BOOST_FOREACH(const string& setname, LHAPDF::availablePDFSets()) {
163
if (!liststr.empty()) liststr += " ";
166
strncpy(s, liststr.c_str(), len);
172
// LHAPDF5 / PDFLIB COMPATIBILITY INTERFACE FUNCTIONS
177
/// LHAPDF library version
178
void getlhapdfversion_(char* s, size_t len) {
179
/// @todo Works? Need to check Fortran string return, string macro treatment, etc.
180
strncpy(s, LHAPDF_VERSION, len);
184
/// Does nothing, only provided for backward compatibility
185
void lhaprint_(int& a) { }
188
/// Set LHAPDF parameters -- does nothing in LHAPDF6!
189
void setlhaparm_(const char* par, int parlength) {
190
/// @todo Can any Fortran LHAPDF params be usefully mapped?
194
/// Return a dummy max number of sets (there is no limitation now)
195
void getmaxnumsets_(int& nmax) {
200
/// Set PDF data path
201
void setpdfpath_(const char* s, size_t len) {
202
/// @todo Works? Need to check C-string copying, null termination
206
LHAPDF::pathsPrepend(s2);
209
/// Get PDF data path (colon-separated if there is more than one element)
210
void getdatapath_(char* s, size_t len) {
211
/// @todo Works? Need to check Fortran string return, string macro treatment, etc.
213
BOOST_FOREACH(const string& path, LHAPDF::paths()) {
214
if (!pathstr.empty()) pathstr += ":";
217
strncpy(s, pathstr.c_str(), len);
221
// PDF initialisation and focus-switching
225
/// @todo Does this version actually take a *path*? What to do?
226
void initpdfsetm_(const int& nset, const char* setpath, int setpathlength) {
227
// Strip file extension for backward compatibility
228
string fullp = string(setpath, setpathlength);
229
// Remove trailing whitespace
230
fullp.erase( std::remove_if( fullp.begin(), fullp.end(), ::isspace ), fullp.end() );
231
// Use only items after the last /
232
const string pap = LHAPDF::dirname(fullp);
233
const string p = LHAPDF::basename(fullp);
234
// Prepend path to search area
235
LHAPDF::pathsPrepend(pap);
237
string path = LHAPDF::file_extn(p).empty() ? p : LHAPDF::file_stem(p);
238
/// @note We correct the misnamed CTEQ6L1/CTEQ6ll set name as a backward compatibility special case.
239
if (boost::algorithm::to_lower_copy(path) == "cteq6ll") path = "cteq6l1";
240
// Create the PDF set with index nset
241
// if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
242
lhapdf_amc::ACTIVESETS[nset] = lhapdf_amc::PDFSetHandler(path); //< @todo Will be wrong if a structured path is given
243
lhapdf_amc::CURRENTSET = nset;
245
/// Load a PDF set (non-multiset version)
246
void initpdfset_(const char* setpath, int setpathlength) {
248
initpdfsetm_(nset1, setpath, setpathlength);
252
/// Load a PDF set by name
253
void initpdfsetbynamem_(const int& nset, const char* setname, int setnamelength) {
254
// Truncate input to size setnamelength
256
p.erase(setnamelength, std::string::npos);
257
// Strip file extension for backward compatibility
258
string name = LHAPDF::file_extn(p).empty() ? p : LHAPDF::file_stem(p);
259
// Remove trailing whitespace
260
name.erase( std::remove_if( name.begin(), name.end(), ::isspace ), name.end() );
261
/// @note We correct the misnamed CTEQ6L1/CTEQ6ll set name as a backward compatibility special case.
262
if (boost::algorithm::to_lower_copy(name) == "cteq6ll") name = "cteq6l1";
263
// Create the PDF set with index nset
264
// if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
265
lhapdf_amc::ACTIVESETS[nset] = lhapdf_amc::PDFSetHandler(name);
266
// Update current set focus
267
lhapdf_amc::CURRENTSET = nset;
269
/// Load a PDF set by name (non-multiset version)
270
void initpdfsetbyname_(const char* setname, int setnamelength) {
272
initpdfsetbynamem_(nset1, setname, setnamelength);
276
/// Load a PDF in current set
277
void initpdfm_(const int& nset, const int& nmember) {
278
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
279
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
280
lhapdf_amc::ACTIVESETS[nset].loadMember(nmember);
281
// Update current set focus
282
lhapdf_amc::CURRENTSET = nset;
284
/// Load a PDF in current set (non-multiset version)
285
void initpdf_(const int& nmember) {
287
initpdfm_(nset1, nmember);
291
/// Get the current set number (i.e. allocation slot index)
292
void getnset_(int& nset) {
293
nset = lhapdf_amc::CURRENTSET;
294
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
295
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
298
/// Explicitly set the current set number (i.e. allocation slot index)
299
void setnset_(const int& nset) {
300
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
301
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
302
lhapdf_amc::CURRENTSET = nset;
306
/// Get the current member number in slot nset
307
void getnmem_(int& nset, int& nmem) {
308
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
309
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
310
nmem = lhapdf_amc::ACTIVESETS[nset].currentmem;
311
// Update current set focus
312
lhapdf_amc::CURRENTSET = nset;
315
/// Set the current member number in slot nset
316
void setnmem_(const int& nset, const int& nmem) {
317
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
318
throw LHAPDF::UserError("Trying to use LHAGLUE set #" +
319
LHAPDF::to_str(nset) + " but it is not initialised");
320
lhapdf_amc::ACTIVESETS[nset].loadMember(nmem);
321
// Update current set focus
322
lhapdf_amc::CURRENTSET = nset;
327
// PDF evolution functions
329
/// Get xf(x) values for common partons from current PDF
330
void evolvepdfm_(const int& nset, const double& x, const double& q, double* fxq) {
331
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
332
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
333
// Evaluate for the 13 LHAPDF5 standard partons (-6..6)
334
for (size_t i = 0; i < 13; ++i) {
336
fxq[i] = lhapdf_amc::ACTIVESETS[nset].activemember()->xfxQ(i-6, x, q);
337
} catch (const exception& e) {
341
// Update current set focus
342
lhapdf_amc::CURRENTSET = nset;
344
/// Get xf(x) values for common partons from current PDF (non-multiset version)
345
void evolvepdf_(const double& x, const double& q, double* fxq) {
347
evolvepdfm_(nset1, x, q, fxq);
350
// PDF evolution functions
351
// NEW BY MZ to evolve one single parton
353
/// Get xf(x) values for common partons from current PDF
354
void evolvepartm_(const int& nset, const int& ipart, const double& x, const double& q, double& fxq) {
355
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
356
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
357
int ipart_copy; // this is to deal with photons, which are labeled 7 in MG5aMC
359
if (ipart==7) ipart_copy = 22;
361
fxq = lhapdf_amc::ACTIVESETS[nset].activemember()->xfxQ(ipart_copy, x, q);
362
} catch (const exception& e) {
365
// Update current set focus
366
lhapdf_amc::CURRENTSET = nset;
368
/// Get xf(x) values for common partons from current PDF (non-multiset version)
369
void evolvepart_( const int& ipart, const double& x, const double& q, double& fxq) {
371
evolvepartm_(nset1, ipart, x, q, fxq);
375
/// Determine if the current PDF has a photon flavour (historically only MRST2004QED)
376
/// @todo Function rather than subroutine?
377
/// @note There is no multiset version. has_photon will respect the current set slot.
379
return lhapdf_amc::ACTIVESETS[lhapdf_amc::CURRENTSET].activemember()->hasFlavor(22);
383
/// Get xfx values from current PDF, including an extra photon flavour
384
void evolvepdfphotonm_(const int& nset, const double& x, const double& q, double* fxq, double& photonfxq) {
385
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
386
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
387
// First evaluate the "normal" partons
388
evolvepdfm_(nset, x, q, fxq);
389
// Then evaluate the photon flavor (historically only for MRST2004QED)
391
photonfxq = lhapdf_amc::ACTIVESETS[nset].activemember()->xfxQ(22, x, q);
392
} catch (const exception& e) {
395
// Update current set focus
396
lhapdf_amc::CURRENTSET = nset;
398
/// Get xfx values from current PDF, including an extra photon flavour (non-multiset version)
399
void evolvepdfphoton_(const double& x, const double& q, double* fxq, double& photonfxq) {
401
evolvepdfphotonm_(nset1, x, q, fxq, photonfxq);
405
/// Get xf(x) values for common partons from a photon PDF
406
void evolvepdfpm_(const int& nset, const double& x, const double& q, const double& p2, const int& ip2, double& fxq) {
407
// Update current set focus
408
lhapdf_amc::CURRENTSET = nset;
409
throw LHAPDF::NotImplementedError("Photon structure functions are not yet supported in LHAPDF6");
411
/// Get xf(x) values for common partons from a photon PDF (non-multiset version)
412
void evolvepdfp_(const double& x, const double& q, const double& p2, const int& ip2, double& fxq) {
414
evolvepdfpm_(nset1, x, q, p2, ip2, fxq);
420
/// Get the alpha_s order for the set
421
void getorderasm_(const int& nset, int& oas) {
422
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
423
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
424
// Set equal to the number of members for the requested set
425
oas = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<int>("AlphaS_OrderQCD");
426
// Update current set focus
427
lhapdf_amc::CURRENTSET = nset;
429
/// Get the alpha_s order for the set (non-multiset version)
430
void getorderas_(int& oas) {
432
getorderasm_(nset1, oas);
436
/// Get the alpha_s(Q) value for set nset
437
double alphaspdfm_(const int& nset, const double& Q){
438
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
439
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
440
return lhapdf_amc::ACTIVESETS[nset].activemember()->alphasQ(Q);
441
// Update current set focus
442
lhapdf_amc::CURRENTSET = nset;
444
/// Get the alpha_s(Q) value for the set (non-multiset version)
445
double alphaspdf_(const double& Q){
447
return alphaspdfm_(nset1, Q);
451
// Metadata functions
453
/// Get the number of error members in the set (with special treatment for single member sets)
454
void numberpdfm_(const int& nset, int& numpdf) {
455
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
456
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
457
// Set equal to the number of members for the requested set
458
numpdf= lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<int>("NumMembers");
459
// Reproduce old LHAPDF v5 behaviour, i.e. subtract 1 if more than 1 member set
460
if (numpdf > 1) numpdf -= 1;
461
// Update current set focus
462
lhapdf_amc::CURRENTSET = nset;
464
/// Get the number of error members in the set (non-multiset version)
465
void numberpdf_(int& numpdf) {
467
numberpdfm_(nset1, numpdf);
471
/// Get the max number of active flavours
472
void getnfm_(const int& nset, int& nf) {
473
//nf = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<int>("AlphaS_NumFlavors");
474
nf = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<int>("NumFlavors");
475
// Update current set focus
476
lhapdf_amc::CURRENTSET = nset;
478
/// Get the max number of active flavours (non-multiset version)
479
void getnf_(int& nf) {
485
/// Get nf'th quark mass
486
void getqmassm_(const int& nset, const int& nf, double& mass) {
487
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
488
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
489
if (nf*nf == 1) mass = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MDown");
490
else if (nf*nf == 4) mass = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MUp");
491
else if (nf*nf == 9) mass = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MStrange");
492
else if (nf*nf == 16) mass = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MCharm");
493
else if (nf*nf == 25) mass = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MBottom");
494
else if (nf*nf == 36) mass = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MTop");
495
else throw LHAPDF::UserError("Trying to get quark mass for invalid quark ID #" + LHAPDF::to_str(nf));
496
// Update current set focus
497
lhapdf_amc::CURRENTSET = nset;
499
/// Get nf'th quark mass (non-multiset version)
500
void getqmass_(const int& nf, double& mass) {
502
getqmassm_(nset1, nf, mass);
506
/// Get the nf'th quark threshold
507
void getthresholdm_(const int& nset, const int& nf, double& Q) {
509
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
510
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
511
if (nf*nf == 1) Q = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdDown");
512
else if (nf*nf == 4) Q = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdUp");
513
else if (nf*nf == 9) Q = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdStrange");
514
else if (nf*nf == 16) Q = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdCharm");
515
else if (nf*nf == 25) Q = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdBottom");
516
else if (nf*nf == 36) Q = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdTop");
517
//else throw LHAPDF::UserError("Trying to get quark threshold for invalid quark ID #" + LHAPDF::to_str(nf));
519
getqmassm_(nset, nf, Q);
521
// Update current set focus
522
lhapdf_amc::CURRENTSET = nset;
524
/// Get the nf'th quark threshold
525
void getthreshold_(const int& nf, double& Q) {
527
getthresholdm_(nset1, nf, Q);
531
/// Print PDF set's description to stdout
532
void getdescm_(const int& nset) {
533
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
534
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
535
cout << lhapdf_amc::ACTIVESETS[nset].activemember()->description() << endl;
536
// Update current set focus
537
lhapdf_amc::CURRENTSET = nset;
545
void getxminm_(const int& nset, const int& nmem, double& xmin) {
546
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
547
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
548
const int activemem = lhapdf_amc::ACTIVESETS[nset].currentmem;
549
lhapdf_amc::ACTIVESETS[nset].loadMember(nmem);
550
xmin = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMin");
551
lhapdf_amc::ACTIVESETS[nset].loadMember(activemem);
552
// Update current set focus
553
lhapdf_amc::CURRENTSET = nset;
555
void getxmin_(const int& nmem, double& xmin) {
557
getxminm_(nset1, nmem, xmin);
561
void getxmaxm_(const int& nset, const int& nmem, double& xmax) {
562
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
563
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
564
const int activemem = lhapdf_amc::ACTIVESETS[nset].currentmem;
565
lhapdf_amc::ACTIVESETS[nset].loadMember(nmem);
566
xmax = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMax");
567
lhapdf_amc::ACTIVESETS[nset].loadMember(activemem);
568
// Update current set focus
569
lhapdf_amc::CURRENTSET = nset;
571
void getxmax_(const int& nmem, double& xmax) {
573
getxmaxm_(nset1, nmem, xmax);
577
void getq2minm_(const int& nset, const int& nmem, double& q2min) {
578
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
579
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
580
const int activemem = lhapdf_amc::ACTIVESETS[nset].currentmem;
581
lhapdf_amc::ACTIVESETS[nset].loadMember(nmem);
582
q2min = LHAPDF::sqr(lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMin"));
583
lhapdf_amc::ACTIVESETS[nset].loadMember(activemem);
584
// Update current set focus
585
lhapdf_amc::CURRENTSET = nset;
587
void getq2min_(const int& nmem, double& q2min) {
589
getq2minm_(nset1, nmem, q2min);
593
void getq2maxm_(const int& nset, const int& nmem, double& q2max) {
594
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
595
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
596
const int activemem = lhapdf_amc::ACTIVESETS[nset].currentmem;
597
lhapdf_amc::ACTIVESETS[nset].loadMember(nmem);
598
q2max = LHAPDF::sqr(lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMax"));
599
lhapdf_amc::ACTIVESETS[nset].loadMember(activemem);
600
// Update current set focus
601
lhapdf_amc::CURRENTSET = nset;
603
void getq2max_(const int& nmem, double& q2max) {
605
getq2maxm_(nset1, nmem, q2max);
609
void getminmaxm_(const int& nset, const int& nmem, double& xmin, double& xmax, double& q2min, double& q2max) {
610
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
611
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
612
const int activemem = lhapdf_amc::ACTIVESETS[nset].currentmem;
613
lhapdf_amc::ACTIVESETS[nset].loadMember(nmem);
614
xmin = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMin");
615
xmax = lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMax");
616
q2min = LHAPDF::sqr(lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMin"));
617
q2max = LHAPDF::sqr(lhapdf_amc::ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMax"));
618
lhapdf_amc::ACTIVESETS[nset].loadMember(activemem);
619
// Update current set focus
620
lhapdf_amc::CURRENTSET = nset;
622
void getminmax_(const int& nmem, double& xmin, double& xmax, double& q2min, double& q2max) {
624
getminmaxm_(nset1, nmem, xmin, xmax, q2min, q2max);
629
/// Backwards compatibility functions for LHAPDF5 calculations of
630
/// PDF uncertainties and PDF correlations (G. Watt, March 2014).
632
// subroutine GetPDFUncTypeM(nset,lMonteCarlo,lSymmetric)
633
void getpdfunctypem_(const int& nset, int& lmontecarlo, int& lsymmetric) {
634
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
635
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
636
const string errorType = lhapdf_amc::ACTIVESETS[nset].activemember()->set().errorType();
637
if (errorType == "replicas") { // Monte Carlo PDF sets
640
} else if (errorType == "symmhessian") { // symmetric eigenvector PDF sets
643
} else { // default: assume asymmetric Hessian eigenvector PDF sets
647
// Update current set focus
648
lhapdf_amc::CURRENTSET = nset;
650
// subroutine GetPDFUncType(lMonteCarlo,lSymmetric)
651
void getpdfunctype_(int& lmontecarlo, int& lsymmetric) {
653
getpdfunctypem_(nset1, lmontecarlo, lsymmetric);
657
// subroutine GetPDFuncertaintyM(nset,values,central,errplus,errminus,errsym)
658
void getpdfuncertaintym_(const int& nset, const double* values, double& central, double& errplus, double& errminus, double& errsymm) {
659
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
660
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
661
const size_t nmem = lhapdf_amc::ACTIVESETS[nset].activemember()->set().size()-1;
662
const vector<double> vecvalues(values, values + nmem + 1);
663
LHAPDF::PDFUncertainty err = lhapdf_amc::ACTIVESETS[nset].activemember()->set().uncertainty(vecvalues, -1);
664
central = err.central;
665
errplus = err.errplus;
666
errminus = err.errminus;
667
errsymm = err.errsymm;
668
// Update current set focus
669
lhapdf_amc::CURRENTSET = nset;
671
// subroutine GetPDFuncertainty(values,central,errplus,errminus,errsym)
672
void getpdfuncertainty_(const double* values, double& central, double& errplus, double& errminus, double& errsymm) {
674
getpdfuncertaintym_(nset1, values, central, errplus, errminus, errsymm);
678
// subroutine GetPDFcorrelationM(nset,valuesA,valuesB,correlation)
679
void getpdfcorrelationm_(const int& nset, const double* valuesA, const double* valuesB, double& correlation) {
680
if (lhapdf_amc::ACTIVESETS.find(nset) == lhapdf_amc::ACTIVESETS.end())
681
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
682
const size_t nmem = lhapdf_amc::ACTIVESETS[nset].activemember()->set().size()-1;
683
const vector<double> vecvaluesA(valuesA, valuesA + nmem + 1);
684
const vector<double> vecvaluesB(valuesB, valuesB + nmem + 1);
685
correlation = lhapdf_amc::ACTIVESETS[nset].activemember()->set().correlation(vecvaluesA,vecvaluesB);
686
// Update current set focus
687
lhapdf_amc::CURRENTSET = nset;
689
// subroutine GetPDFcorrelation(valuesA,valuesB,correlation)
690
void getpdfcorrelation_(const double* valuesA, const double* valuesB, double& correlation) {
692
getpdfcorrelationm_(nset1, valuesA, valuesB, correlation);
696
///////////////////////////////////////
699
/// REALLY OLD PDFLIB COMPATILITY FUNCTIONS
701
/// PDFLIB initialisation function
702
void pdfset_(const char* par, const double* value, int parlength) {
704
// Identify the calling program (yuck!)
706
if (my_par.find("NPTYPE") != string::npos) {
707
cout << "==== LHAPDF6 USING PYTHIA-TYPE LHAGLUE INTERFACE ====" << endl;
708
// Take PDF ID from value[2]
709
lhapdf_amc::ACTIVESETS[1] = lhapdf_amc::PDFSetHandler(value[2]+1000*value[1]);
710
} else if (my_par.find("HWLHAPDF") != string::npos) {
711
cout << "==== LHAPDF6 USING HERWIG-TYPE LHAGLUE INTERFACE ====" << endl;
712
// Take PDF ID from value[0]
713
lhapdf_amc::ACTIVESETS[1] = lhapdf_amc::PDFSetHandler(value[0]);
714
} else if (my_par.find("DEFAULT") != string::npos) {
715
cout << "==== LHAPDF6 USING DEFAULT-TYPE LHAGLUE INTERFACE ====" << endl;
716
// Take PDF ID from value[0]
717
lhapdf_amc::ACTIVESETS[1] = lhapdf_amc::PDFSetHandler(value[0]);
719
cout << "==== LHAPDF6 USING PDFLIB-TYPE LHAGLUE INTERFACE ====" << endl;
720
// Take PDF ID from value[2]
721
lhapdf_amc::ACTIVESETS[1] = lhapdf_amc::PDFSetHandler(value[2]+1000*value[1]);
724
lhapdf_amc::CURRENTSET = 1;
726
// Extract parameters for common blocks (with sensible fallback values)
727
lhapdf_amc::PDFPtr pdf = lhapdf_amc::ACTIVESETS[1].activemember();
728
w50513_.xmin = pdf->info().get_entry_as<double>("XMin", 0.0);
729
w50513_.xmax = pdf->info().get_entry_as<double>("XMax", 1.0);
730
w50513_.q2min = LHAPDF::sqr(pdf->info().get_entry_as<double>("QMin", 1.0));
731
w50513_.q2max = LHAPDF::sqr(pdf->info().get_entry_as<double>("QMax", 1.0e5));
732
w50512_.qcdl4 = pdf->info().get_entry_as<double>("AlphaS_Lambda4", 0.0);
733
w50512_.qcdl5 = pdf->info().get_entry_as<double>("AlphaS_Lambda5", 0.0);
734
lhapdfr_.qcdlha4 = pdf->info().get_entry_as<double>("AlphaS_Lambda4", 0.0);
735
lhapdfr_.qcdlha5 = pdf->info().get_entry_as<double>("AlphaS_Lambda5", 0.0);
737
// Activate legacy/compatibility LHAPDF5-type behaviour re. broken Lambda values
738
if (pdf->info().get_entry_as<bool>("Pythia6LambdaV5Compat", true)) {
739
w50512_.qcdl4 = 0.192;
740
w50512_.qcdl5 = 0.192;
741
lhapdfr_.qcdlha4 = 0.192;
742
lhapdfr_.qcdlha5 = 0.192;
746
/// PDFLIB nucleon structure function querying
747
void structm_(const double& x, const double& q,
748
double& upv, double& dnv, double& usea, double& dsea,
749
double& str, double& chm, double& bot, double& top, double& glu) {
750
lhapdf_amc::CURRENTSET = 1;
751
/// Fill (partial) parton return variables
752
lhapdf_amc::PDFPtr pdf = lhapdf_amc::ACTIVESETS[1].activemember();
753
dsea = pdf->xfxQ(-1, x, q);
754
usea = pdf->xfxQ(-2, x, q);
755
dnv = pdf->xfxQ(1, x, q) - dsea;
756
upv = pdf->xfxQ(2, x, q) - usea;
757
str = pdf->xfxQ(3, x, q);
758
chm = (pdf->hasFlavor(4)) ? pdf->xfxQ(4, x, q) : 0;
759
bot = (pdf->hasFlavor(5)) ? pdf->xfxQ(5, x, q) : 0;
760
top = (pdf->hasFlavor(6)) ? pdf->xfxQ(6, x, q) : 0;
761
glu = pdf->xfxQ(21, x, q);
764
/// PDFLIB photon structure function querying
765
void structp_(const double& x, const double& q2, const double& p2, const double& ip2,
766
double& upv, double& dnv, double& usea, double& dsea,
767
double& str, double& chm, double& bot, double& top, double& glu) {
768
throw LHAPDF::NotImplementedError("Photon structure functions are not yet supported");
771
/// PDFLIB statistics on PDF under/overflows
773
/// @note Can't do anything...
780
// LHAPDF namespace C++ compatibility code
781
#ifdef ENABLE_LHAGLUE_CXX
784
void LHAPDF::setVerbosity(LHAPDF::Verbosity noiselevel) {
785
LHAPDF::setVerbosity((int) noiselevel);
788
void LHAPDF::setPDFPath(const string& path) {
792
string LHAPDF::pdfsetsPath() {
796
int LHAPDF::numberPDF() {
801
int LHAPDF::numberPDF(int nset) {
803
numberpdfm_(nset,nmem);
807
void LHAPDF::initPDF( int memset) {
809
initpdfm_(nset1, memset);
811
void LHAPDF::initPDF(int nset, int memset) {
812
initpdfm_(nset, memset);
816
double LHAPDF::xfx(double x, double Q, int fl) {
817
vector<double> r(13);
818
evolvepdf_(x, Q, &r[0]);
821
double LHAPDF::xfx(int nset, double x, double Q, int fl) {
822
vector<double> r(13);
823
evolvepdfm_(nset, x, Q, &r[0]);
827
vector<double> LHAPDF::xfx(double x, double Q) {
828
vector<double> r(13);
829
evolvepdf_(x, Q, &r[0]);
832
vector<double> LHAPDF::xfx(int nset, double x, double Q) {
833
vector<double> r(13);
834
evolvepdfm_(nset, x, Q, &r[0]);
838
void LHAPDF::xfx(double x, double Q, double* results) {
839
evolvepdf_(x, Q, results);
841
void LHAPDF::xfx(int nset, double x, double Q, double* results) {
842
evolvepdfm_(nset, x, Q, results);
846
vector<double> LHAPDF::xfxphoton(double x, double Q) {
847
vector<double> r(13);
849
evolvepdfphoton_(x, Q, &r[0], mphoton);
850
r.push_back(mphoton);
853
vector<double> LHAPDF::xfxphoton(int nset, double x, double Q) {
854
vector<double> r(13);
856
evolvepdfphotonm_(nset, x, Q, &r[0], mphoton);
857
r.push_back(mphoton);
861
void LHAPDF::xfxphoton(double x, double Q, double* results) {
862
evolvepdfphoton_(x, Q, results, results[13]);
864
void LHAPDF::xfxphoton(int nset, double x, double Q, double* results) {
865
evolvepdfphotonm_(nset, x, Q, results, results[13]);
868
double LHAPDF::xfxphoton(double x, double Q, int fl) {
869
vector<double> r(13);
871
evolvepdfphoton_(x, Q, &r[0], mphoton);
872
if (fl == 7) return mphoton;
875
double LHAPDF::xfxphoton(int nset, double x, double Q, int fl) {
876
vector<double> r(13);
878
evolvepdfphotonm_(nset, x, Q, &r[0], mphoton);
879
if ( fl == 7 ) return mphoton;
884
void LHAPDF::initPDFSet(const string& filename, int nmem) {
885
initPDFSet(1,filename, nmem);
888
void LHAPDF::initPDFSet(int nset, const string& filename, int nmem) {
889
initPDFSetByName(nset,filename);
890
ACTIVESETS[nset].loadMember(nmem);
895
void LHAPDF::initPDFSet(const string& filename, SetType type ,int nmem) {
896
// silently ignore type
897
initPDFSet(1,filename, nmem);
900
void LHAPDF::initPDFSet(int nset, const string& filename, SetType type ,int nmem) {
901
// silently ignore type
902
initPDFSetByName(nset,filename);
903
ACTIVESETS[nset].loadMember(nmem);
907
void LHAPDF::initPDFSet(int nset, int setid, int nmem) {
908
ACTIVESETS[nset] = PDFSetHandler(setid); //
912
void LHAPDF::initPDFSet(int setid, int nmem) {
913
initPDFSet(1,setid,nmem);
917
void LHAPDF::initPDFSetByName(const string& filename) {
918
std::cout << "initPDFSetByName: " << filename << std::endl;
919
char cfilename[SIZE+1];
920
strncpy(cfilename, filename.c_str(), SIZE);
921
initpdfsetbyname_(cfilename, filename.length());
924
void LHAPDF::initPDFSetByName(int nset, const string& filename) {
925
char cfilename[SIZE+1];
926
strncpy(cfilename, filename.c_str(), SIZE);
927
initpdfsetbynamem_(nset, cfilename, filename.length());
930
void LHAPDF::initPDFSetByName(const string& filename, SetType type) {
931
//silently ignore type
932
std::cout << "initPDFSetByName: " << filename << std::endl;
933
char cfilename[SIZE+1];
934
strncpy(cfilename, filename.c_str(), SIZE);
935
initpdfsetbyname_(cfilename, filename.length());
938
void LHAPDF::initPDFSetByName(int nset, const string& filename, SetType type) {
939
//silently ignore type
940
char cfilename[SIZE+1];
941
strncpy(cfilename, filename.c_str(), SIZE);
942
initpdfsetbynamem_(nset, cfilename, filename.length());
946
void LHAPDF::getDescription() {
950
void LHAPDF::getDescription(int nset) {
951
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
952
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
953
cout << ACTIVESETS[nset].activemember()->set().description() << endl;
957
double LHAPDF::alphasPDF(double Q) {
958
return LHAPDF::alphasPDF(1, Q) ;
961
double LHAPDF::alphasPDF(int nset, double Q) {
962
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
963
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
965
// return alphaS for the requested set
966
return ACTIVESETS[nset].activemember()->alphasQ(Q);
970
bool LHAPDF::hasPhoton(){
971
return has_photon_();
975
int LHAPDF::getOrderAlphaS() {
976
return LHAPDF::getOrderAlphaS(1) ;
979
int LHAPDF::getOrderAlphaS(int nset) {
980
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
981
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
983
// return alphaS Order for the requested set
984
return ACTIVESETS[nset].activemember()->info().get_entry_as<int>("AlphaS_OrderQCD", -1);
988
int LHAPDF::getOrderPDF() {
989
return LHAPDF::getOrderPDF(1) ;
992
int LHAPDF::getOrderPDF(int nset) {
993
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
994
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
996
// return PDF order for the requested set
997
return ACTIVESETS[nset].activemember()->info().get_entry_as<int>("OrderQCD", -1);
1001
double LHAPDF::getLam4(int nmem) {
1002
return LHAPDF::getLam4(1, nmem) ;
1005
double LHAPDF::getLam4(int nset, int nmem) {
1006
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
1007
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
1009
ACTIVESETS[nset].loadMember(nmem);
1010
return ACTIVESETS[nset].activemember()->info().get_entry_as<double>("AlphaS_Lambda4", -1.0);
1014
double LHAPDF::getLam5(int nmem) {
1015
return LHAPDF::getLam5(1, nmem) ;
1018
double LHAPDF::getLam5(int nset, int nmem) {
1019
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
1020
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
1022
ACTIVESETS[nset].loadMember(nmem);
1023
return ACTIVESETS[nset].activemember()->info().get_entry_as<double>("AlphaS_Lambda5", -1.0);
1027
int LHAPDF::getNf() {
1028
return LHAPDF::getNf(1) ;
1031
int LHAPDF::getNf(int nset) {
1032
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
1033
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
1035
// return alphaS Order for the requested set
1036
return ACTIVESETS[nset].activemember()->info().get_entry_as<int>("NumFlavors");
1040
double LHAPDF::getXmin(int nmem) {
1041
return LHAPDF::getXmin(1, nmem) ;
1044
double LHAPDF::getXmin(int nset, int nmem) {
1045
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
1046
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
1048
// return alphaS Order for the requested set
1049
ACTIVESETS[nset].loadMember(nmem);
1050
return ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMin");
1053
double LHAPDF::getXmax(int nmem) {
1054
return LHAPDF::getXmax(1, nmem) ;
1057
double LHAPDF::getXmax(int nset, int nmem) {
1058
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
1059
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
1061
// return alphaS Order for the requested set
1062
ACTIVESETS[nset].loadMember(nmem);
1063
return ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMax");
1066
double LHAPDF::getQ2min(int nmem) {
1067
return LHAPDF::getQ2min(1, nmem) ;
1070
double LHAPDF::getQ2min(int nset, int nmem) {
1071
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
1072
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
1074
// return alphaS Order for the requested set
1075
ACTIVESETS[nset].loadMember(nmem);
1076
return pow(ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMin"),2);
1079
double LHAPDF::getQ2max(int nmem) {
1080
return LHAPDF::getQ2max(1,nmem) ;
1083
double LHAPDF::getQ2max(int nset, int nmem) {
1084
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
1085
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
1087
// return alphaS Order for the requested set
1088
ACTIVESETS[nset].loadMember(nmem);
1089
return pow(ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMax"),2);
1092
double LHAPDF::getQMass(int nf) {
1093
return LHAPDF::getQMass(1, nf) ;
1096
double LHAPDF::getQMass(int nset, int nf) {
1098
getqmassm_(nset, nf, mass);
1102
double LHAPDF::getThreshold(int nf) {
1103
return LHAPDF::getThreshold(1, nf) ;
1106
double LHAPDF::getThreshold(int nset, int nf) {
1108
getthresholdm_(nset, nf, thres);
1112
void LHAPDF::usePDFMember(int member) {
1116
void LHAPDF::usePDFMember(int nset, int member) {
1117
initpdfm_(nset, member);
1120
#endif // ENABLE_LHAGLUE_CXX