diff --git a/src/pdf.C b/src/pdf.C index 4d3df8a..78276e5 100644 --- a/src/pdf.C +++ b/src/pdf.C @@ -1,229 +1,273 @@ #include "pdf.h" #include "dyres_interface.h" #include "interface.h" #include "scales.h" #include "settings.h" #include #include #include LHAPDF::PDF* pdf::lhapdf = 0; +void (*pdf::xfxq)(const double &x, const double &Q, double *fPDF) = 0; +double (*pdf::alphas)(const double &Q) = 0; + +void pdf::lhaxfxq(const double &x, const double &Q, double *r) +{ + vector fPDF; fPDF.resize(13); + lhapdf->xfxQ(x,Q,fPDF); + copy(fPDF.begin(),fPDF.end(),r); +} + +double pdf::lhaalphas(const double &Q) +{ + return lhapdf->alphasQ(Q); +} + void pdf::init() { - //printf(" ==Initialize PDF set from LHAPDF==\n\n"); - //printf("\n"); - LHAPDF::Info& cfg = LHAPDF::getConfig(); - - cfg.set_entry("Verbosity" , 0 ); - cfg.set_entry("Interpolator" , "logcubic" ); - cfg.set_entry("Extrapolator" , "continuation" ); - cfg.set_entry("ForcePositive" , 0 ); - cfg.set_entry("AlphaS_Type" , "analytic" ); - cfg.set_entry("MZ" , 91.1876 ); - cfg.set_entry("MUp" , 0.002 ); - cfg.set_entry("MDown" , 0.005 ); - cfg.set_entry("MStrange" , 0.10 ); - cfg.set_entry("MCharm" , 1.29 ); - cfg.set_entry("MBottom" , 4.19 ); - cfg.set_entry("MTop" , 172.9 ); - cfg.set_entry("Pythia6LambdaV5Compat" , true ); - - //Old interface - LHAPDF::initPDFSet(opts.LHAPDFset); //LHAPDF::initPDFSetByName(opts.LHAPDFset); - LHAPDF::initPDF(opts.LHAPDFmember); - - //New interface - lhapdf = LHAPDF::mkPDF(opts.LHAPDFset, opts.LHAPDFmember); - - if (opts.PDFerrors && LHAPDF::numberPDF() > 1) + if (!opts.externalpdf) { - opts.totpdf = LHAPDF::numberPDF()+1; - pdferropts_.pdferr_ = true; - pdferropts_.totpdf_ = LHAPDF::numberPDF()+1; + //printf(" ==Initialize PDF set from LHAPDF==\n\n"); + //printf("\n"); + LHAPDF::Info& cfg = LHAPDF::getConfig(); + + cfg.set_entry("Verbosity" , 0 ); + cfg.set_entry("Interpolator" , "logcubic" ); + cfg.set_entry("Extrapolator" , "continuation" ); + cfg.set_entry("ForcePositive" , 0 ); + cfg.set_entry("AlphaS_Type" , "analytic" ); + cfg.set_entry("MZ" , 91.1876 ); + cfg.set_entry("MUp" , 0.002 ); + cfg.set_entry("MDown" , 0.005 ); + cfg.set_entry("MStrange" , 0.10 ); + cfg.set_entry("MCharm" , 1.29 ); + cfg.set_entry("MBottom" , 4.19 ); + cfg.set_entry("MTop" , 172.9 ); + cfg.set_entry("Pythia6LambdaV5Compat" , true ); + + //Old interface + LHAPDF::initPDFSet(opts.LHAPDFset); //LHAPDF::initPDFSetByName(opts.LHAPDFset); + LHAPDF::initPDF(opts.LHAPDFmember); + + //New interface + lhapdf = LHAPDF::mkPDF(opts.LHAPDFset, opts.LHAPDFmember); + + xfxq = lhaxfxq; + alphas = lhaalphas; + + if (opts.PDFerrors && LHAPDF::numberPDF() > 1) + { + opts.totpdf = LHAPDF::numberPDF()+1; + pdferropts_.pdferr_ = true; + pdferropts_.totpdf_ = LHAPDF::numberPDF()+1; + } + else + { + opts.totpdf = 1; + pdferropts_.pdferr_ = false; + pdferropts_.totpdf_ = 1; + } } else { opts.totpdf = 1; pdferropts_.pdferr_ = false; pdferropts_.totpdf_ = 1; } // initialization of alphas setalphas(); // read g from the PDF setg(); //take the cmass and b mass from the PDF // cmass=dsqrt(mcsq) // bmass=dsqrt(mbsq) } //set value of alphas void pdf::setalphas() { + //Old LHAPDF interface //couple_.amz_=LHAPDF::alphasPDF(dymasses_.zmass_); - couple_.amz_ = lhapdf->alphasQ(dymasses_.zmass_); + + //New LHAPDF interface + //couple_.amz_ = lhapdf->alphasQ(dymasses_.zmass_); + + //Allows external alphas + couple_.amz_ = alphas(dymasses_.zmass_); + double scale = fabs(scale_.scale_); if (opts_.approxpdf_ == 1) { int nloop = 3; qcdcouple_.as_=dyalphas_mcfm_(scale,couple_.amz_,nloop); } else //qcdcouple_.as_=dyalphas_lhapdf_(scale); - //qcdcouple_.as_=LHAPDF::alphasPDF(scale); - qcdcouple_.as_=lhapdf->alphasQ(scale); + //qcdcouple_.as_=LHAPDF::alphasPDF(scale); //Old LHAPDF interface + //qcdcouple_.as_=lhapdf->alphasQ(scale); //New LHAPDF interface + qcdcouple_.as_=alphas(scale); //Allows external alphas qcdcouple_.ason2pi_=qcdcouple_.as_/(2*M_PI); qcdcouple_.ason4pi_=qcdcouple_.as_/(4*M_PI); qcdcouple_.gsq_=4*M_PI*qcdcouple_.as_; } //set the value of the g-parameter of the non perturbative form factor void pdf::setg() { LHAPDF::PDFInfo info(opts.LHAPDFset, opts.LHAPDFmember); double gformfactor = info.get_entry_as("g", -1); if (gformfactor >= 0) { cout << "g form factor: input from PDF member: " << gformfactor << endl; opts.g_param = gformfactor; g_param_.g_param_ = gformfactor; np_.g_ = opts.g_param; } } void dysetpdf_(int& member) { + if (opts.externalpdf) return; + if (member == 0) { if (opts.PDFerrors && opts.totpdf > 1) LHAPDF::initPDF(0); else LHAPDF::initPDF(opts.LHAPDFmember); } else LHAPDF::initPDF(member); pdf::setalphas(); // setg(); //set g separately when setting a different PDF in the resummed part } void setmellinpdf_(int &member){ // if member is still same than dont do anything //if (lastMember==member) return; //if (v_mellinpdf.size() 1. || x <= 0.) { for (int i = -MAXNF; i <= MAXNF; i++) fx[MAXNF+i]=0.; return; } - + + //Old LHAPDF interface //double fPDF[13]; //LHAPDF::xfx(x,xmu,fPDF); - vector fPDF; fPDF.resize(13); - pdf::lhapdf->xfxQ(x,xmu,fPDF); + //New LHAPDF interface + //vector fPDF; fPDF.resize(13); + //pdf::lhapdf->xfxQ(x,xmu,fPDF); + + //Allows external PDF + double fPDF[13]; + pdf::xfxq(x,xmu,fPDF); + //vector pids = pdf::lhapdf->flavors(); //for (int i = 0; i < pids.size(); i++) //cout << pids[i] << endl; // cout << endl; // for (int i = -MAXNF; i <= MAXNF; i++) // cout << fPDF[6+i]/x << " "; // cout << endl; //xmu *= xmu; //fPDF[6-5] = pdf::lhapdf->xfxQ(-5,x,xmu); //fPDF[6-4] = pdf::lhapdf->xfxQ(-4,x,xmu); //fPDF[6-3] = pdf::lhapdf->xfxQ(-3,x,xmu); //fPDF[6-2] = pdf::lhapdf->xfxQ(-2,x,xmu); //fPDF[6-1] = pdf::lhapdf->xfxQ(-1,x,xmu); //fPDF[6+0] = pdf::lhapdf->xfxQ(21,x,xmu); //fPDF[6+1] = pdf::lhapdf->xfxQ(1 ,x,xmu); //fPDF[6+2] = pdf::lhapdf->xfxQ(2 ,x,xmu); //fPDF[6+3] = pdf::lhapdf->xfxQ(3 ,x,xmu); //fPDF[6+4] = pdf::lhapdf->xfxQ(4 ,x,xmu); //fPDF[6+5] = pdf::lhapdf->xfxQ(5 ,x,xmu); // for (int i = -MAXNF; i <= MAXNF; i++) // cout << fPDF[6+i]/x << " "; // cout << endl; if (ih == 1) //proton for (int i = -MAXNF; i <= MAXNF; i++) fx[MAXNF+i]=fPDF[6+i]/x; else if (ih == -1) //antiproton for (int i = -MAXNF; i <= MAXNF; i++) fx[MAXNF+i]=fPDF[6-i]/x; //switch off flavours //fx[MAXNF-5]=0.; //bbar //fx[MAXNF-4]=0.; //cbar //fx[MAXNF-3]=0.; //sbar //fx[MAXNF-2]=0.; //fx[MAXNF-1]=0.; //fx[MAXNF+0]=0.; //gluon //fx[MAXNF+1]=0.; //fx[MAXNF+2]=0.; //fx[MAXNF+3]=0.; //s //fx[MAXNF+4]=0.; //c //fx[MAXNF+5]=0.; //b //impose positivity //fx[MAXNF-5]=max(0.,fx[MAXNF-5]); //fx[MAXNF-4]=max(0.,fx[MAXNF-4]); //fx[MAXNF-3]=max(0.,fx[MAXNF-3]); //fx[MAXNF-2]=max(0.,fx[MAXNF-2]); //fx[MAXNF-1]=max(0.,fx[MAXNF-1]); //fx[MAXNF+0]=max(0.,fx[MAXNF+0]); //fx[MAXNF+1]=max(0.,fx[MAXNF+1]); //fx[MAXNF+2]=max(0.,fx[MAXNF+2]); //fx[MAXNF+3]=max(0.,fx[MAXNF+3]); //fx[MAXNF+4]=max(0.,fx[MAXNF+4]); //fx[MAXNF+5]=max(0.,fx[MAXNF+5]); //make u = d and ubar = dbar //fx[MAXNF-1]=fx[MAXNF-2]; //fx[MAXNF+1]=fx[MAXNF+2]; /* if (x < 90./13000.*exp(-1.)) { fx[MAXNF-5]=0.; fx[MAXNF-4]=0.; fx[MAXNF-3]=0.; fx[MAXNF-2]=0.; fx[MAXNF-1]=0.; fx[MAXNF+0]=0.; fx[MAXNF+1]=0.; fx[MAXNF+2]=0.; fx[MAXNF+3]=0.; fx[MAXNF+4]=0.; fx[MAXNF+5]=0.; } */ } diff --git a/src/pdf.h b/src/pdf.h index 6a643a5..c77d407 100644 --- a/src/pdf.h +++ b/src/pdf.h @@ -1,25 +1,31 @@ #ifndef pdf_h #define pdf_h #include #include "mcfm_interface.h" namespace pdf { extern LHAPDF::PDF* lhapdf; extern void init(); extern void setalphas(); extern void setg(); + + extern void (*xfxq)(const double &x, const double &Q, double *fPDF); + extern void lhaxfxq(const double &x, const double &Q, double *fPDF); + + extern double (*alphas)(const double &Q); + extern double lhaalphas(const double &Q); } extern "C" { void fdist_(int& ih, double& x, double& xmu, double fx[2*MAXNF+1]); void dysetpdf_(int& member); void setmellinpdf_(int& member); } #endif