diff --git a/src/pdf.C b/src/pdf.C index 78276e5..dfa711d 100644 --- a/src/pdf.C +++ b/src/pdf.C @@ -1,273 +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() { if (!opts.externalpdf) { //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_); //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); //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; + opts.g1 = gformfactor; g_param_.g_param_ = gformfactor; - np_.g_ = opts.g_param; + np_.g_ = opts.g1; } } 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); //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.; } */ }