From 844125d7b29120cc93d15ad7a042d92f1d47df67 Mon Sep 17 00:00:00 2001 From: Toni Makela <toni.makela@cern.ch> Date: Tue, 31 Aug 2021 15:29:45 +0200 Subject: [PATCH 1/3] Added CIJET reaction --- CMakeLists.txt | 1 + reactions/CIJET/CMakeLists.txt | 11 + reactions/CIJET/include/CIJETReader.h | 73 +++ reactions/CIJET/include/ReactionCIJET.h | 47 ++ reactions/CIJET/src/CIJETReader.cc | 105 ++++ reactions/CIJET/src/ReactionCIJET.cc | 135 +++++ reactions/CIJET/src/cijet.f | 554 +++++++++++++++++++++ reactions/CIJET/yaml/parameters.yaml | 6 + reactions/CIJET/yaml/parameters_fixed.yaml | 6 + reactions/CIJET/yaml/parameters_orig.yaml | 8 + 10 files changed, 946 insertions(+) create mode 100644 reactions/CIJET/CMakeLists.txt create mode 100644 reactions/CIJET/include/CIJETReader.h create mode 100644 reactions/CIJET/include/ReactionCIJET.h create mode 100644 reactions/CIJET/src/CIJETReader.cc create mode 100644 reactions/CIJET/src/ReactionCIJET.cc create mode 100644 reactions/CIJET/src/cijet.f create mode 100644 reactions/CIJET/yaml/parameters.yaml create mode 100644 reactions/CIJET/yaml/parameters_fixed.yaml create mode 100644 reactions/CIJET/yaml/parameters_orig.yaml diff --git a/CMakeLists.txt b/CMakeLists.txt index e0f57e645..b05eb5052 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -148,6 +148,7 @@ add_subdirectory(reactions/HVQMNR_LHCb_7TeV_charm) add_subdirectory(reactions/Hathor) add_subdirectory(reactions/HathorSingleTop) add_subdirectory(reactions/fastNLO) +add_subdirectory(reactions/CIJET) #Tools add_subdirectory(tools/draw/src) diff --git a/reactions/CIJET/CMakeLists.txt b/reactions/CIJET/CMakeLists.txt new file mode 100644 index 000000000..5df9e3be2 --- /dev/null +++ b/reactions/CIJET/CMakeLists.txt @@ -0,0 +1,11 @@ +set(TARGET reactionCIJET) +add_library(${TARGET} MODULE + src/ReactionCIJET.cc + src/CIJETReader.cc + src/cijet.f +) +target_link_libraries(${TARGET} PRIVATE xfitter) +target_include_directories(${TARGET} PUBLIC include) +install(TARGETS ${TARGET} DESTINATION ${DEST_MODULE}) +install(FILES yaml/parameters.yaml DESTINATION ${DEST_YAML}/reactions RENAME ReactionCIJET.yaml) + diff --git a/reactions/CIJET/include/CIJETReader.h b/reactions/CIJET/include/CIJETReader.h new file mode 100644 index 000000000..3b3af46ba --- /dev/null +++ b/reactions/CIJET/include/CIJETReader.h @@ -0,0 +1,73 @@ +#ifndef CIJET_H +#define CIJET_H + +//-------------------------------------------------------------- +// C++ wraper for using CIJET interface + + +#include <vector> +#include <string> +#include <iomanip> +#include <iostream> +#include <fstream> +#include <math.h> +#include <cstdlib> +#include "ReactionTheory.h" +//-------------------------------------------- + +using namespace std; + +class CIJETReader { + + public: + +// initialization + CIJETReader(string fname) { + filepds=fname; + } + + CIJETReader(string fname, ReactionTheory* reaction) { + filepds=fname; + _reactionTheory=reaction; + } + +// name of table file + string filepds; + +// QCD parameters + int Oqcd=1; + double mur=1.0, muf=1.0; + +// contact interaction parameters + int cshape=1; + double lambda, c1, c2, c3; + +// normalization + double fut=1.; + +// setting parameters + void setnorm(double norm); + void setorder(int od); + void setscales(double mr, double mf); + void setci(double lam, vector<double> coe); + +// initialization + void setinit(); + +// calculating xsecs + vector<double> calcxsec(); + + private: + +// commons + ReactionTheory* _reactionTheory; + +// Get strong coupling + double GetAs_(double mr); + +// Get PDFs (recalculate, e.g. call to xFitter routine) + void GetPDF_(double x, double mf, double pdf[13]); + +}; + +#endif diff --git a/reactions/CIJET/include/ReactionCIJET.h b/reactions/CIJET/include/ReactionCIJET.h new file mode 100644 index 000000000..a70396827 --- /dev/null +++ b/reactions/CIJET/include/ReactionCIJET.h @@ -0,0 +1,47 @@ +#ifndef xFitter_ReactionCIJET +#define xFitter_ReactionCIJET + +#pragma once + +#include "ReactionTheory.h" +#include "CIJETReader.h" +//#include <memory> +#include <map> +#include <vector> +#include "hf_errlog.h" + +/** + @class' ReactionCIJET + + @brief A wrapper class for CIJET reaction + + Based on the ReactionTheory class + + @version 0.3 + @date 2016-12-06 + */ + +class CIJETReaction : public CIJETReader { +public: + CIJETReaction(string name, ReactionTheory* reaction) : CIJETReader(name, reaction) {}; +protected: + CIJETReaction(string name) : CIJETReader(name) {}; // not public! +}; + +class ReactionCIJET : public ReactionTheory { +public: + ReactionCIJET(){}; + + +public: + virtual string getReactionName() const { return "CIJET" ;}; + virtual void initTerm(TermData* td) override final; + virtual void atStart() {} //< nothing todo + virtual void compute(TermData*, valarray<double> &val, map<string, valarray<double> > &err); +protected: + virtual int parseOptions(){ return 0;}; + // CIJET inherited functions + std::map<int,std::vector<CIJETReaction*> > ffnlo; +}; + +#endif diff --git a/reactions/CIJET/src/CIJETReader.cc b/reactions/CIJET/src/CIJETReader.cc new file mode 100644 index 000000000..a25895018 --- /dev/null +++ b/reactions/CIJET/src/CIJETReader.cc @@ -0,0 +1,105 @@ +#include <string> +#include <cstring> +#include "CIJETReader.h" +#include "ReactionTheory.h" + +using namespace std; + +// Interface to xFitter FORTRAN routines +extern "C" +{ + void cijetinit_(const char* fname, double* mufip, double* murip, int* od, double* norm); + void cijetxsec_(const char* fname, double* lam, double* cl, int* npt, double* pres); + // wrappers for applgrid functionality + double appl_fnalphas_(double* Q); + void appl_fnpdf_(double* x, double* mf, double* pdf); +} +double appl_fnalphas_(double* Q) +{ + return alphas_wrapper_(*Q); +} +void appl_fnpdf_(double* x, double* mf, double* pdf) +{ + double q = *mf; + std::valarray<double> pdfV(13); + pdf_xfxq_wrapper_(*x, q, &pdfV[0]); + for(int i = 0; i < 13; i++) pdf[i] = pdfV[i]; +} + +//------------------------------------------------------------------------------------ + +void CIJETReader::setinit() +{ + //note the character lenth in Fortran is hardwired to 100 + char ff[101]; + strncpy(ff, filepds.c_str(), 100); + cijetinit_(ff, &muf, &mur, &Oqcd, &fut); +} + +vector<double> CIJETReader::calcxsec(void) +{ + char ff[101]; + strncpy(ff, filepds.c_str(), 100); + double invlamsq, cpl[3]; + double pres[200]; + int npt; + + invlamsq=1.0/lambda/lambda; + cpl[0]=c1; + cpl[1]=c2; + cpl[2]=c3; + + cijetxsec_(ff, &invlamsq, cpl, &npt, pres); + + vector<double> xsec; + for (int i=0; i<npt; i++) { +// cout<<"gaojun "<<pres[i]<<endl; + xsec.push_back(pres[i]); + } + return xsec; +} + +void CIJETReader::setnorm(double norm) +{ + fut=norm; // overall normalization factor +// cout<<"=========> CInorm set to "<<fut<<endl; +} + +void CIJETReader::setorder(int od) +{ + Oqcd=od; //0 or 1 for LO or NLO +// cout<<"=========> CIorder set to "<<Oqcd<<endl; +} + +void CIJETReader::setscales(double mf, double mr) +{ + muf=mf; + mur=mr; +// cout<<"=========> CIscales set to "<<muf<<" "<<mur<<endl; +} + +void CIJETReader::setci(double lam, vector<double> coe) +{ + lambda=lam; // Lambda in TeV + c1=coe[0]; // color singlet coefficient, LL + c2=coe[1]; //LR + c3=coe[2]; //RR +} + +void CIJETReader::GetPDF_(double x, double mf, double pdf[13]) +{ + double q = mf; + std::valarray<double> pdfV(13); + //_reactionTheory->xfx(x, q, &pdfV[0]); //Not supported in xFitter 2.1 + pdf_xfxq_wrapper_(x, q, &pdfV[0]); + for(int i = 0; i < 13; i++) + pdf[i] = pdfV[i]; +} + +double CIJETReader::GetAs_(double mr) +{ + double q = mr; + //return _reactionTheory->alphaS(q); //Not supported in xFitter 2.1 + return alphas_wrapper_(q); +} + diff --git a/reactions/CIJET/src/ReactionCIJET.cc b/reactions/CIJET/src/ReactionCIJET.cc new file mode 100644 index 000000000..49ae13c4c --- /dev/null +++ b/reactions/CIJET/src/ReactionCIJET.cc @@ -0,0 +1,135 @@ +// DB 08/2017 +/* + @file ReactionCIJET.cc + @date 2016-12-06 + @author AddReaction.py + Created by AddReaction.py on 2016-12-06 +*/ + +#include "ReactionCIJET.h" + + +//______________________________________________________________________________ +// the class factories +extern "C" ReactionCIJET* create() { + return new ReactionCIJET(); +} +//______________________________________________________________________________ +// Initialise for a given dataset: +void ReactionCIJET::initTerm(TermData* td) { + + int ID = td->id; // ID=dataSetID -> updated to termdata td + if ( td->hasParam("Filename") ) { + // --- Read file and instantiate CIJET + const std::string filename = td->getParamS("Filename"); + hf_errlog(17086501,"I: Reading CIJET file "+filename); + ffnlo.insert(std::make_pair(ID,vector<CIJETReaction*>{new CIJETReaction(filename,this)} )); + } + else if ( td->hasParam("Filenames") ) { + const std::string filenames = td->getParamS("Filenames"); + std::stringstream ss(filenames); + std::string filename; + while (std::getline(ss, filename, ',')) { + hf_errlog(17086501,"I: Reading CIJET file "+filename); + ffnlo[ID].push_back(new CIJETReaction(filename,this)); + } + } + else { + hf_errlog(17086503,"F:No CIJET file specified. Please provide 'Filename' or 'Filenames' to reaction CIJET."); + return; + } + + for ( CIJETReaction* fnlo : ffnlo[ID] ) { + // --- Set order of calculation + if ( td->hasParam("Order") ) { // Local order + hf_errlog(17094510,"W: Ignoring key 'Order' in .dat file. Only global parameter 'Order' is used."); + } + string order = td->getParamS("Order"); // Global order + bool success=true; + // CIJET default is 'NLO' + if (order=="NNLO" ) fnlo->setorder(1); // swith NLO ON + else if (order=="NLO" ) fnlo->setorder(1); // swith NLO ON + else if (order=="LO" ) fnlo->setorder(0); // switch NLO OFF + else hf_errlog(17094502,"E: CIJET. Unrecognized order: "+order); + if (!success) hf_errlog(17094503,"W: CIJET. Requested order "+order+" cannot be set."); + + // --- Set CI normalization + double nm=1.0; + if ( td->hasParam("CInorm") ) { // Local order + nm = *td->getParamD("CInorm"); + } + fnlo->setnorm(nm); + + // --- Set scale factors + double cmur=1, cmuf=1; + if ( td->hasParam("ScaleFacMuR") ) { + hf_errlog(17094504,"I: Setting CIJET scale factor mu_R: "+td->getParamS("ScaleFacMuR")); + cmur = *td->getParamD("ScaleFacMuR"); + } + //else { + // cmur=GetParam("ScaleFacMuR"); + //} + if ( td->hasParam("ScaleFacMuF") ) { + hf_errlog(17094505,"I: Setting CIJET scale factor mu_F: "+td->getParamS("ScaleFacMuF")); + cmuf = *td->getParamD("ScaleFacMuF"); + } + //else { + // cmuf=GetParam("ScaleFacMuF"); + //} + fnlo->setscales(cmuf,cmur); + + // ---Read of the table + fnlo->setinit(); + + } +} + +//______________________________________________________________________________ +// Main function to compute results at an iteration +void ReactionCIJET::compute(TermData* td, valarray<double> &val, map<string, valarray<double> > &err) +{ + td->actualizeWrappers(); + // Get relevant parameters: + int cicase = td->getParamI("CIcase"); + double lamotev = *td->getParamD("LamoTeV"); + double cp1 = *td->getParamD("CI1"); + double cp3 = *td->getParamD("CI3"); + double cp5 = *td->getParamD("CI5"); + // Adjust according to type of contact interactions + // 1-3 constrained fit + if ( cicase==1 ) { //pure left or right handed + cp3=0.0; //using cp1 as free parameter + cp5=0.0; + } + else if ( cicase==2 ) { //vector case + cp3=2*cp1; //using cp1 as free parameter + cp5=cp1; + } + else if ( cicase==3 ) { //axial-vector case + cp3=-2*cp1; //using cp1 as free parameter + cp5=cp1; + } + else { //default free-N parameters + } + + vector<double> ci; + ci.push_back(cp1); + ci.push_back(cp3); + ci.push_back(cp5); + + vector<double> valfnlo; + valfnlo.reserve(val.size()); + int ID = td->id; + for ( CIJETReaction* fnlo : ffnlo[ID] ) { + fnlo->setci(lamotev,ci); + std::vector<double> cs = fnlo->calcxsec(); + for ( std::size_t i=0; i<cs.size(); i++) valfnlo.push_back(cs[i]); + } + + if ( val.size()!=valfnlo.size() ) + hf_errlog(17116501,"F: Size of CIJET cross section array does not match data."); + for ( std::size_t i=0; i<val.size(); i++) val[i]=(i<valfnlo.size())?valfnlo[i]:0.; +} + + + diff --git a/reactions/CIJET/src/cijet.f b/reactions/CIJET/src/cijet.f new file mode 100644 index 000000000..cbe8f50ae --- /dev/null +++ b/reactions/CIJET/src/cijet.f @@ -0,0 +1,554 @@ +c CIJET1.1 interface made for Xfitter, Jun Gao, 2018.09.10 + +c performing grid reading and initialization + subroutine cijetinit(fname, mufip, murip, oqcd, fut) + +! new line + implicit none +! maximum number of data sets and points per set + integer, parameter :: mset=1, mpint=200 + integer, parameter :: xnd=40, qnd=15, nbk=42 + +! input scale choices, CI shape choice, only color-singlet + real(8) :: mufip, murip, fut + integer :: cshape ! 1 for LL, 2 for RR, 3 for vector, 4 for av + real(8) :: cpl(3)=0.d0 ! color singlet couplings, LL, LR, RR + real(8) :: ash(10) ! individual coupling combinations +! input QCD order + integer :: oqcd ! 0 for LO, 1 for NLO +! grid name + character*100 :: fname + + integer, save :: befirst=1 +! storing total number of sets + integer :: cset + +! storing reference Lambda scale, no of points, fname per set + real(8) :: Nscale(mset) ! hardwired currently + integer :: pset(mset) + character*100 :: ffnm(mset) +! storing pp or ppbar, and mur/muf + integer :: sig(mset) + real(8) :: srof(mset) + +! storing reference x-Q-mu scale per point + real(8) :: xgd(mset,mpint,xnd) + real(8) :: qgd(mset,mpint,qnd) + real(8) :: smu(mset,mpint) + integer :: nxgd(mset,mpint) + integer :: nqgd(mset,mpint) + +! storing grid weight per point + real(8) :: ciwgt(10,xnd,xnd*qnd,5,mset,mpint) + +! beta function needed for scale evaluation + real(8), parameter :: beta0=(11*3.d0-2*5.d0)/12.d0/3.1415926d0 + real(8), parameter :: scaf(9)=(/0.5,0.5,0.5,1.,1.,1.,2.,2.,2./) + real(8), parameter :: scar(9)=(/0.5,1.,2.,0.5,1.,2.,0.5,1.,2./) + +! individual set workspace + integer :: nt, cnt, pointor(3*nbk*mpint+1) + integer :: i, j, k, lhc, ang, scalescheme, Nx, Nq + real(8) :: chid, chiu, masscutd, masscutu, mu0, Rcone, fct + integer :: sid, bid1, bid2, nstart, nend, nca + real(8) :: sqrts, slog, wgt(5), as, sc, tm1, tm2, tm3, tm4 + character*100 :: xx + +! shared common block for latter convolution + common /cigrid/ Nscale, xgd, qgd, smu, ciwgt, + - srof, cset, pset, nxgd, nqgd, sig, ffnm + +! initialization of commons + if(befirst==1) then + sid=0;cset=0;pset=0;srof=1.d0 + ciwgt=0.d0;Nscale=5000.d0 + xgd=0.d0;qgd=0.d0;smu=0.d0 + befirst=0 + endif + + write(*,*) "*****************************************" + write(*,*) "* Reading CIJET1.1 table with parameters:" + write(*,*) "* mufip, murip, oqcd, fut, fname=" + write(*,*) "*", mufip, murip, oqcd, fut, trim(fname) + write(*,*) "* Jun Gao, 2018.09.10" + write(*,*) "*****************************************" + +! end of initialization of commons + +! define the couplings apart from 1/Lambda^2 +! select case (cshape) +! case(3) +! cpl(1)=1.d0;cpl(2)= 2.d0;cpl(3)=1.d0 +! case(4) +! cpl(1)=1.d0;cpl(2)=-2.d0;cpl(3)=1.d0 +! case default +! cpl(1)=1.d0;cpl(2)=0.d0;cpl(3)=0.d0 +! end select + ash(1) =1.d0 !cpl(1)+cpl(3) + ash(2) =1.d0 !cpl(1)^2+cpl(3)^2 + ash(3) =1.d0 !cpl(2)^2 + ash(4) =1.d0 !cpl(1)+cpl(3) + ash(5) =1.d0 !cpl(2) + ash(6) =1.d0 !cpl(1)+cpl(3) + ash(7) =1.d0 !cpl(2) + ash(8) =1.d0 !cpl(1)^2+cpl(3)^2 + ash(9) =1.d0 !cpl(2)^2 + ash(10)=1.d0 !cpl(1)^2+cpl(3)^2 + +! first load to find index of components + nt=1 + cnt=0 + pointor=0 + open(unit=155, file=fname, status="old", action="read") + do while(0.lt.1) + read(155, *, end=200) xx + cnt=cnt+1 + if(xx.eq."JCC") then + pointor(nt)=cnt + nt=nt+1 + endif + enddo + 200 pointor(nt)=cnt+1 + + if(mod(nt-1, 3*nbk).ne.0) then + print *, "wrong grid file read!" + stop + endif + close(155) + +! determine the correct scale choice + do i=1, 9 + if(mufip==scaf(i)) then + sid=1+(i-1)/3 + endif + enddo + if(sid==0) then + print *, "scale not available!";stop + endif + slog=dlog(murip) + +! update record + cset=cset+1 + pset(cset)=(nt-1)/(3*nbk) + ffnm(cset)=fname + sig(cset)=1 + srof(cset)=murip/mufip +! end first + + open(unit=156, file=fname, status="old", action="read") +! begin of the main loop on the kinematic bins + do nca=1, pset(cset) + + read(156, *) lhc, xx, sc, as + read(156, *) sqrts, ang, scalescheme, Rcone + read(156, *) chid, chiu, masscutd, masscutu, mu0 + read(156, *) Nx, Nq + +! transfer bin-xsec to double differential + fct=fut/(chiu-chid)/(masscutu-masscutd) + +! storing needed grid information + if((Nx.gt.xnd).or.(Nq.gt.qnd)) then + print *, "X or Q grid nodes wrong!"; stop + endif + + nxgd(cset,nca)=Nx + nqgd(cset,nca)=Nq + + do i=1, Nx + read(156, *) xgd(cset,nca,i) + enddo + do i=1, Nq + read(156, *) qgd(cset,nca,i) + enddo + smu(cset,nca)=mu0 + if(lhc==0) sig(cset)=-1 + + +! must include all color-singlet case + if(xx.eq."fitll") then + print *, "grid file should be for CS at least!" + stop + endif + +! loop over scale choices + do j=1, 3 + +! loop over 42 coefficients + do k=1, nbk + +! skip to the correct scale choice + bid1=nbk*(3*(nca-1)+j-1)+k + bid2=nbk*(3*(nca-1)+j-1)+k+1 + + nstart=pointor(bid1) + nend=pointor(bid2)-1 + +! account for the header lines + if((j==3).and.(k==nbk).and.(nca.ne.pset(cset))) then + nend=nend-(pointor(1)-1) + endif + +! drop JCC + read(156, *) xx + +! block by block + do i=nstart+1, nend +! skipping + if(j.ne.sid) then + read(156, *) xx + endif +! working + if(j.eq.sid) then + read(156, *) nx, nq, wgt(1:5) +! tranfer to double differentail + wgt=wgt*fct +! adding all components with correct couplings + select case (k) +! O(0) +! inteference + case (1) +! ciwgt(1,nx,nq,1:5,cset,nca)=ciwgt(1,nx,nq,1:5,cset,nca) +! - +(cpl(1)+cpl(3))*wgt(1:5) + ciwgt(1,nx,nq,1:5,cset,nca)=ciwgt(1,nx,nq,1:5,cset,nca) + - +(ash(1))*wgt(1:5) +! square + case (5) +! ciwgt(2,nx,nq,1:5,cset,nca)=ciwgt(2,nx,nq,1:5,cset,nca) +! - +(cpl(1)**2+cpl(3)**2)*wgt(1:5) + ciwgt(2,nx,nq,1:5,cset,nca)=ciwgt(2,nx,nq,1:5,cset,nca) + - +(ash(2))*wgt(1:5) + case (8) +! ciwgt(2,nx,nq,1:5,cset,nca)=ciwgt(2,nx,nq,1:5,cset,nca) +! - +(cpl(2)**2)*wgt(1:5) + ciwgt(3,nx,nq,1:5,cset,nca)=ciwgt(3,nx,nq,1:5,cset,nca) + - +(ash(3))*wgt(1:5) +! O(1) +! inteference, const + case (11) + if(oqcd.ne.0) then +! ciwgt(3,nx,nq,1:5,cset,nca)=ciwgt(3,nx,nq,1:5,cset,nca) +! - +(cpl(1)+cpl(3))*wgt(1:5) + ciwgt(4,nx,nq,1:5,cset,nca)=ciwgt(4,nx,nq,1:5,cset,nca) + - +(ash(4))*wgt(1:5) + endif + case (15) + if(oqcd.ne.0) then +! ciwgt(3,nx,nq,1:5,cset,nca)=ciwgt(3,nx,nq,1:5,cset,nca) +! - +(cpl(2))*wgt(1:5) + ciwgt(5,nx,nq,1:5,cset,nca)=ciwgt(5,nx,nq,1:5,cset,nca) + - +(ash(5))*wgt(1:5) + endif +! inteference, logs + case (12) + if(oqcd.ne.0) then +! ciwgt(4,nx,nq,1:5,cset,nca)=ciwgt(4,nx,nq,1:5,cset,nca) +! - +(cpl(1)+cpl(3))*wgt(1:5) + ciwgt(6,nx,nq,1:5,cset,nca)=ciwgt(6,nx,nq,1:5,cset,nca) + - +(ash(6))*wgt(1:5) + endif + case (16) + if(oqcd.ne.0) then +! ciwgt(4,nx,nq,1:5,cset,nca)=ciwgt(4,nx,nq,1:5,cset,nca) +! - +(cpl(2))*wgt(1:5) + ciwgt(7,nx,nq,1:5,cset,nca)=ciwgt(7,nx,nq,1:5,cset,nca) + - +(ash(7))*wgt(1:5) + endif +! square, const + case (19) + if(oqcd.ne.0) then +! ciwgt(5,nx,nq,1:5,cset,nca)=ciwgt(5,nx,nq,1:5,cset,nca) +! - +(cpl(1)**2+cpl(3)**2)*wgt(1:5) + ciwgt(8,nx,nq,1:5,cset,nca)=ciwgt(8,nx,nq,1:5,cset,nca) + - +(ash(8))*wgt(1:5) + endif + case (25) + if(oqcd.ne.0) then +! ciwgt(5,nx,nq,1:5,cset,nca)=ciwgt(5,nx,nq,1:5,cset,nca) +! - +(cpl(2)**2)*wgt(1:5) + ciwgt(9,nx,nq,1:5,cset,nca)=ciwgt(9,nx,nq,1:5,cset,nca) + - +(ash(9))*wgt(1:5) + endif +! square, logs + case (20) + if(oqcd.ne.0) then +! ciwgt(6,nx,nq,1:5,cset,nca)=ciwgt(6,nx,nq,1:5,cset,nca) +! - +(cpl(1)**2+cpl(3)**2)*wgt(1:5) + ciwgt(10,nx,nq,1:5,cset,nca)=ciwgt(10,nx,nq,1:5,cset,nca) + - +(ash(10))*wgt(1:5) + endif + case default + read(156, *) xx + end select + + endif + + enddo +! end of block + + enddo +! end of coefficients + +! adding back mur dependence + if((j.eq.sid).and.(oqcd.ne.0)) then + ciwgt(4,1:xnd,1:qnd*xnd,1:5,cset,nca)= + - ciwgt(4,1:xnd,1:qnd*xnd,1:5,cset,nca) + - +2.d0*slog*beta0*ciwgt(1,1:xnd,1:qnd*xnd,1:5,cset,nca) + endif + + enddo +! end of scale choices + + enddo +! end of kinematic loop + + close(156) + + return + + end subroutine cijetinit + +c perform grid convolution and calculation of xsecs + subroutine cijetxsec(fname, invlamsq, cpl, npt, pres) + +! new line + implicit none +! maximum number of data sets and points per set + integer, parameter :: mset=1, mpint=200 + integer, parameter :: xnd=40, qnd=15, nbk=42 +! input 1/Lambda^2 [in TeV^-2] + real(8) :: invlamsq, res(mpint), pres(mpint) +! sign (+ for destructive, -for con.) + real(8) :: cpl(3) ! color singlet couplings, LL, LR, RR + real(8) :: ash(10) ! individual coupling combinations +! number of data points and additional multi-factor + integer :: npt +! grid name + character*100 :: fname + +! storing total number of sets + integer :: cset + +! storing reference Lambda scale, no of points, fname per set + real(8) :: Nscale(mset) ! hardwired currently + integer :: pset(mset) + character*100 :: ffnm(mset) +! storing pp or ppbar, and mur/muf + integer :: sig(mset) + real(8) :: srof(mset) + +! storing reference x-Q-mu scale per point + real(8) :: xgd(mset,mpint,xnd) + real(8) :: qgd(mset,mpint,qnd) + real(8) :: smu(mset,mpint) + integer :: nxgd(mset,mpint) + integer :: nqgd(mset,mpint) + +! storing grid weight per point + real(8) :: ciwgt(10,xnd,xnd*qnd,5,mset,mpint) + +! workspace + real(8) :: XPDF(-6:6, xnd, qnd), xi(xnd), qi(qnd) + real(8) :: pdfpg(xnd*qnd, xnd, 5) + integer :: pid, nca, Nx, Nq, sgg + integer :: i, j, k, ii, jj, kk, iqd + real(8) :: norm(qnd), tmp, cialphas, tm1, tm2, tm3, tm4 + external :: cialphas, cievolPDF + + +! shared common block for latter convolution + common /cigrid/ Nscale, xgd, qgd, smu, ciwgt, + - srof, cset, pset, nxgd, nqgd, sig, ffnm + +! find the correct ID + pid=0; res=0.d0 + do i=1, cset + if(fname==ffnm(i)) pid=i + enddo + if(pid==0) then + print *, "CI grid not loaded!",fname; stop + endif + +!test print *, "gaojun1 !", invlamsq, cpl, cset, pid, pset(pid) + +! find the correct couplings + ash(1) =cpl(1)+cpl(3) + ash(2) =cpl(1)**2+cpl(3)**2 + ash(3) =cpl(2)**2 + ash(4) =cpl(1)+cpl(3) + ash(5) =cpl(2) + ash(6) =cpl(1)+cpl(3) + ash(7) =cpl(2) + ash(8) =cpl(1)**2+cpl(3)**2 + ash(9) =cpl(2)**2 + ash(10)=cpl(1)**2+cpl(3)**2 + +! loop over the kinematic bins + do nca=1, pset(pid) + + XPDF=0.d0 +! passing grid configuration + Nx=nxgd(pid,nca); Nq=nqgd(pid,nca); sgg=sig(pid) + xi=xgd(pid, nca, 1:xnd); qi=qgd(pid, nca, 1:qnd) + +! calculate the PDFs in each grid points + do i=1, Nx + do j=1, Nq + call cievolPDF(xi(i), qi(j), XPDF(-6, i, j)) + do k=-6, 6 +! interpolation kernal is x**2*pdf + XPDF(k, i, j)=XPDF(k, i, j)*xi(i) + enddo + enddo + enddo +! end of PDF calls + +! calculation of the PDF grid quantity + pdfpg=0.d0 + + do i=1, Nx + do ii=1, Nx + do j=1, Nq + do k=1, 4 + do kk=k+1, 5 + pdfpg(j+Nq*(ii-1), i, 1)= + - pdfpg(j+Nq*(ii-1), i, 1) + - +XPDF(k, i, j)*XPDF(sgg*kk, ii, j) + - +XPDF(kk, i, j)*XPDF(sgg*k, ii, j) + - +XPDF(-k, i, j)*XPDF(-sgg*kk, ii, j) + - +XPDF(-kk, i, j)*XPDF(-sgg*k, ii, j) + pdfpg(j+Nq*(ii-1), i, 3)= + - pdfpg(j+Nq*(ii-1), i, 3) + - +XPDF(k, i, j)*XPDF(-sgg*kk, ii, j) + - +XPDF(kk, i, j)*XPDF(-sgg*k, ii, j) + - +XPDF(sgg*k, i, j)*XPDF(-kk, ii, j) + - +XPDF(sgg*kk, i, j)*XPDF(-k, ii, j) + enddo + enddo + do k=1, 5 + pdfpg(j+Nq*(ii-1), i, 2)= + - pdfpg(j+Nq*(ii-1), i, 2) + - +XPDF(k, i, j)*XPDF(sgg*k, ii, j) + - +XPDF(-k, i, j)*XPDF(-sgg*k, ii, j) + pdfpg(j+Nq*(ii-1), i, 4)= + - pdfpg(j+Nq*(ii-1), i, 4) + - +XPDF(k, i, j)*XPDF(-sgg*k, ii, j) + - +XPDF(sgg*k, i, j)*XPDF(-k, ii, j) + pdfpg(j+Nq*(ii-1), i, 5)= + - pdfpg(j+Nq*(ii-1), i, 5) + - +XPDF(k, i, j)*XPDF(0, ii, j) + - +XPDF(sgg*k, i, j)*XPDF(0, ii, j) + - +XPDF(-k, i, j)*XPDF(0, ii, j) + - +XPDF(-sgg*k, i, j)*XPDF(0, ii, j) + enddo + enddo + enddo + enddo +! end PDF grid + +! calculate alphas grid quantity + do i=1, Nq + norm(i)=cialphas(srof(pid)*qi(i)) + enddo +! end alphas grid + +! now ready for cross sections +! inteference + tmp=0.d0 + do k=1, 5 + do ii=1, Nx + do j=1, Nq + do i=1, Nx + iqd=j+(ii-1)*Nq +! O(0) + tmp=tmp+pdfpg(iqd,i,k)*(ash(1)*ciwgt(1,i,iqd,k,pid,nca)) + - *norm(j) +! O(1) const + tmp=tmp+pdfpg(iqd,i,k)*(ash(4)*ciwgt(4,i,iqd,k,pid,nca)+ + - ash(5)*ciwgt(5,i,iqd,k,pid,nca)) + - *norm(j)**2 +! O(1) logs + tmp=tmp+pdfpg(iqd,i,k)*(ash(6)*ciwgt(6,i,iqd,k,pid,nca)+ + - ash(7)*ciwgt(7,i,iqd,k,pid,nca)) + - *(-0.5d0*dlog(abs(invlamsq)*(smu(pid,nca)/1.0d3)**2)) + - *norm(j)**2 + enddo + enddo + enddo + enddo + res(nca)=res(nca)+tmp*((Nscale(pid)/1.d3)**2)*invlamsq +! square + tmp=0.d0 + do k=1, 5 + do ii=1, Nx + do j=1, Nq + do i=1, Nx + iqd=j+(ii-1)*Nq +! O(0) + tmp=tmp+pdfpg(iqd,i,k)*(ash(2)*ciwgt(2,i,iqd,k,pid,nca)+ + - ash(3)*ciwgt(3,i,iqd,k,pid,nca)) +! O(1) const + tmp=tmp+pdfpg(iqd,i,k)*(ash(8)*ciwgt(8,i,iqd,k,pid,nca)+ + - ash(9)*ciwgt(9,i,iqd,k,pid,nca)) + - *norm(j) +! O(1) logs + tmp=tmp+pdfpg(iqd,i,k)*(ash(10)*ciwgt(10,i,iqd,k,pid,nca)) + - *(-0.5d0*dlog(abs(invlamsq)*(smu(pid,nca)/1.0d3)**2)) + - *norm(j) + enddo + enddo + enddo + enddo + res(nca)=res(nca)+tmp*(((Nscale(pid)/1.d3)**2)*invlamsq)**2 +! end of cross sections + + enddo +! end of loop on kinematic bins + + npt=pset(pid); pres=0.d0 + pres(1:pset(pid))=res(1:pset(pid)) + + return + +c end of PDF calculation + end subroutine cijetxsec + + +!!!!!!must rely on external Fortran routines for alphas +!!!!!!and PDFs; can not really wrapped them into reaction +!!!!!!class due to problem of passing functions within class + +c link to PDF routines + subroutine cievolPDF(x, q, res) + +! new line + implicit none + real(8) :: x, q, res(-6:6) + real(8) :: xf(-6:6) + + res=0.d0 +! e.g., + call appl_fnpdf(x, Q, xf) +! test with LHAPDF call evolvePDF(x, Q, xf) + + res(-6:6)=xf + return + + end subroutine cievolPDF + +c link to PDF routines + real(8) function cialphas(q) + +! new line + implicit none + real(8) :: q + real(8), external :: appl_fnalphas +! test with LHAPDF real(8), external :: alphasPDF + + cialphas=0.118d0 +! e.g., + cialphas=appl_fnalphas(Q) +! test with LHAPDF cialphas=alphasPDF(Q) + return + + end function cialphas diff --git a/reactions/CIJET/yaml/parameters.yaml b/reactions/CIJET/yaml/parameters.yaml new file mode 100644 index 000000000..c707be4be --- /dev/null +++ b/reactions/CIJET/yaml/parameters.yaml @@ -0,0 +1,6 @@ +CIJET : + LamoTeV : 10.0 + CIcase : 1 +CI1 : -0.50 +CI3 : 0.00 +CI5 : 0.00 diff --git a/reactions/CIJET/yaml/parameters_fixed.yaml b/reactions/CIJET/yaml/parameters_fixed.yaml new file mode 100644 index 000000000..b4e990246 --- /dev/null +++ b/reactions/CIJET/yaml/parameters_fixed.yaml @@ -0,0 +1,6 @@ +CIJET : + LamoTeV : 9.0 + CIcase : 1 +CI1 : 0.00 +CI3 : 0.00 +CI5 : 0.00 diff --git a/reactions/CIJET/yaml/parameters_orig.yaml b/reactions/CIJET/yaml/parameters_orig.yaml new file mode 100644 index 000000000..0b297fe3a --- /dev/null +++ b/reactions/CIJET/yaml/parameters_orig.yaml @@ -0,0 +1,8 @@ +CIJET : + LamoTeV : 10.0 + CIcase : 1 +CI1 : + value: 1.00 + step: 0.01 +CI3 : 1.00 +CI5 : 1.00 -- GitLab From 4a7696beedc4453f3a84ea23faf62a885fedb8f8 Mon Sep 17 00:00:00 2001 From: Toni Makela <toni.makela@cern.ch> Date: Thu, 18 Nov 2021 14:10:27 +0100 Subject: [PATCH 2/3] Reset to original status, added up-to-date CIJET reaction files, now about to checkout, pull & merge master to CIJET before pushing to git --- reactions/CIJET/include/CIJETReader.h | 51 +++-- reactions/CIJET/include/ReactionCIJET.h | 18 +- reactions/CIJET/src/CIJETReader.cc | 96 +++------- reactions/CIJET/src/ReactionCIJET.cc | 209 ++++++++++----------- reactions/CIJET/src/cijet.f | 25 +-- reactions/CIJET/yaml/parameters_fixed.yaml | 6 - reactions/CIJET/yaml/parameters_orig.yaml | 8 - 7 files changed, 169 insertions(+), 244 deletions(-) delete mode 100644 reactions/CIJET/yaml/parameters_fixed.yaml delete mode 100644 reactions/CIJET/yaml/parameters_orig.yaml diff --git a/reactions/CIJET/include/CIJETReader.h b/reactions/CIJET/include/CIJETReader.h index 3b3af46ba..b71fa717f 100644 --- a/reactions/CIJET/include/CIJETReader.h +++ b/reactions/CIJET/include/CIJETReader.h @@ -13,60 +13,55 @@ #include <math.h> #include <cstdlib> #include "ReactionTheory.h" +#include "TermData.h" + //-------------------------------------------- using namespace std; class CIJETReader { - public: +public: // initialization - CIJETReader(string fname) { - filepds=fname; - } - - CIJETReader(string fname, ReactionTheory* reaction) { - filepds=fname; + CIJETReader(string fname) { + filepds=fname; + } + + CIJETReader(string fname, ReactionTheory* reaction) { + filepds=fname; _reactionTheory=reaction; - } + } // name of table file - string filepds; + string filepds; // QCD parameters - int Oqcd=1; - double mur=1.0, muf=1.0; + int Oqcd=1; + double mur=1.0, muf=1.0; // contact interaction parameters - int cshape=1; - double lambda, c1, c2, c3; + double lambda, c1, c2, c3; // normalization - double fut=1.; + double fut=1.; // setting parameters - void setnorm(double norm); - void setorder(int od); - void setscales(double mr, double mf); - void setci(double lam, vector<double> coe); + void setnorm(double norm); + void setorder(int od); + void setscales(double mr, double mf); + void setci(double lam, vector<double> coe); // initialization - void setinit(); + void setinit(); // calculating xsecs - vector<double> calcxsec(); + vector<double> calcxsec(); - private: +private: // commons - ReactionTheory* _reactionTheory; - -// Get strong coupling - double GetAs_(double mr); - -// Get PDFs (recalculate, e.g. call to xFitter routine) - void GetPDF_(double x, double mf, double pdf[13]); + ReactionTheory* _reactionTheory; }; diff --git a/reactions/CIJET/include/ReactionCIJET.h b/reactions/CIJET/include/ReactionCIJET.h index a70396827..80c762541 100644 --- a/reactions/CIJET/include/ReactionCIJET.h +++ b/reactions/CIJET/include/ReactionCIJET.h @@ -23,9 +23,9 @@ class CIJETReaction : public CIJETReader { public: - CIJETReaction(string name, ReactionTheory* reaction) : CIJETReader(name, reaction) {}; + CIJETReaction(string name, ReactionTheory* reaction) : CIJETReader(name, reaction) {}; protected: - CIJETReaction(string name) : CIJETReader(name) {}; // not public! + CIJETReaction(string name) : CIJETReader(name) {}; // not public! }; class ReactionCIJET : public ReactionTheory { @@ -34,14 +34,14 @@ public: public: - virtual string getReactionName() const { return "CIJET" ;}; - virtual void initTerm(TermData* td) override final; - virtual void atStart() {} //< nothing todo - virtual void compute(TermData*, valarray<double> &val, map<string, valarray<double> > &err); + virtual string getReactionName() const { return "CIJET" ;}; + virtual void initTerm(TermData* td) override final; + virtual void atStart() {} //< nothing todo + virtual void compute(TermData*, valarray<double> &val, map<string, valarray<double> > &err); protected: - virtual int parseOptions(){ return 0;}; - // CIJET inherited functions - std::map<int,std::vector<CIJETReaction*> > ffnlo; + virtual int parseOptions(){ return 0;}; + // CIJET inherited functions + std::map<int,std::vector<CIJETReaction*> > ffnlo; }; #endif diff --git a/reactions/CIJET/src/CIJETReader.cc b/reactions/CIJET/src/CIJETReader.cc index a25895018..57eb691fa 100644 --- a/reactions/CIJET/src/CIJETReader.cc +++ b/reactions/CIJET/src/CIJETReader.cc @@ -8,98 +8,64 @@ using namespace std; // Interface to xFitter FORTRAN routines extern "C" { - void cijetinit_(const char* fname, double* mufip, double* murip, int* od, double* norm); - void cijetxsec_(const char* fname, double* lam, double* cl, int* npt, double* pres); - // wrappers for applgrid functionality - double appl_fnalphas_(double* Q); - void appl_fnpdf_(double* x, double* mf, double* pdf); -} -double appl_fnalphas_(double* Q) -{ - return alphas_wrapper_(*Q); -} -void appl_fnpdf_(double* x, double* mf, double* pdf) -{ - double q = *mf; - std::valarray<double> pdfV(13); - pdf_xfxq_wrapper_(*x, q, &pdfV[0]); - for(int i = 0; i < 13; i++) pdf[i] = pdfV[i]; + void cijetinit_(const char* fname, double* mufip, double* murip, int* od, double* norm); + void cijetxsec_(const char* fname, double* lam, double* cl, int* npt, double* pres); } //------------------------------------------------------------------------------------ void CIJETReader::setinit() { - //note the character lenth in Fortran is hardwired to 100 - char ff[101]; - strncpy(ff, filepds.c_str(), 100); - cijetinit_(ff, &muf, &mur, &Oqcd, &fut); + //note the character length in Fortran is hardwired to 100 + char ff[101]; + strncpy(ff, filepds.c_str(), 100); + cijetinit_(ff, &muf, &mur, &Oqcd, &fut); } vector<double> CIJETReader::calcxsec(void) { - char ff[101]; - strncpy(ff, filepds.c_str(), 100); - double invlamsq, cpl[3]; - double pres[200]; - int npt; + char ff[101]; + strncpy(ff, filepds.c_str(), 100); + double invlamsq, cpl[3]; + double pres[200]; + int npt; + + invlamsq=1.0/lambda/lambda; + cpl[0]=c1; + cpl[1]=c2; + cpl[2]=c3; - invlamsq=1.0/lambda/lambda; - cpl[0]=c1; - cpl[1]=c2; - cpl[2]=c3; + cijetxsec_(ff, &invlamsq, cpl, &npt, pres); - cijetxsec_(ff, &invlamsq, cpl, &npt, pres); + vector<double> xsec; + for (int i=0; i<npt; i++) xsec.push_back(pres[i]); - vector<double> xsec; - for (int i=0; i<npt; i++) { -// cout<<"gaojun "<<pres[i]<<endl; - xsec.push_back(pres[i]); - } - return xsec; + return xsec; } void CIJETReader::setnorm(double norm) { - fut=norm; // overall normalization factor -// cout<<"=========> CInorm set to "<<fut<<endl; + fut=norm; // overall normalization factor +// cout<<"=========> CInorm set to "<<fut<<endl; } void CIJETReader::setorder(int od) { - Oqcd=od; //0 or 1 for LO or NLO -// cout<<"=========> CIorder set to "<<Oqcd<<endl; + Oqcd=od; //0 or 1 for LO or NLO +// cout<<"=========> CIorder set to "<<Oqcd<<endl; } void CIJETReader::setscales(double mf, double mr) { - muf=mf; - mur=mr; -// cout<<"=========> CIscales set to "<<muf<<" "<<mur<<endl; + muf=mf; + mur=mr; +// cout<<"=========> CIscales set to "<<muf<<" "<<mur<<endl; } void CIJETReader::setci(double lam, vector<double> coe) { - lambda=lam; // Lambda in TeV - c1=coe[0]; // color singlet coefficient, LL - c2=coe[1]; //LR - c3=coe[2]; //RR -} - -void CIJETReader::GetPDF_(double x, double mf, double pdf[13]) -{ - double q = mf; - std::valarray<double> pdfV(13); - //_reactionTheory->xfx(x, q, &pdfV[0]); //Not supported in xFitter 2.1 - pdf_xfxq_wrapper_(x, q, &pdfV[0]); - for(int i = 0; i < 13; i++) - pdf[i] = pdfV[i]; -} - -double CIJETReader::GetAs_(double mr) -{ - double q = mr; - //return _reactionTheory->alphaS(q); //Not supported in xFitter 2.1 - return alphas_wrapper_(q); + lambda=lam; // Lambda in TeV + c1=coe[0]; // color singlet coefficient, LL + c2=coe[1]; //LR + c3=coe[2]; //RR } - diff --git a/reactions/CIJET/src/ReactionCIJET.cc b/reactions/CIJET/src/ReactionCIJET.cc index 49ae13c4c..8fc0e0398 100644 --- a/reactions/CIJET/src/ReactionCIJET.cc +++ b/reactions/CIJET/src/ReactionCIJET.cc @@ -12,124 +12,113 @@ //______________________________________________________________________________ // the class factories extern "C" ReactionCIJET* create() { - return new ReactionCIJET(); + return new ReactionCIJET(); } //______________________________________________________________________________ // Initialise for a given dataset: void ReactionCIJET::initTerm(TermData* td) { - - int ID = td->id; // ID=dataSetID -> updated to termdata td - if ( td->hasParam("Filename") ) { - // --- Read file and instantiate CIJET - const std::string filename = td->getParamS("Filename"); - hf_errlog(17086501,"I: Reading CIJET file "+filename); - ffnlo.insert(std::make_pair(ID,vector<CIJETReaction*>{new CIJETReaction(filename,this)} )); - } - else if ( td->hasParam("Filenames") ) { - const std::string filenames = td->getParamS("Filenames"); - std::stringstream ss(filenames); - std::string filename; - while (std::getline(ss, filename, ',')) { - hf_errlog(17086501,"I: Reading CIJET file "+filename); - ffnlo[ID].push_back(new CIJETReaction(filename,this)); - } - } - else { - hf_errlog(17086503,"F:No CIJET file specified. Please provide 'Filename' or 'Filenames' to reaction CIJET."); - return; - } - - for ( CIJETReaction* fnlo : ffnlo[ID] ) { - // --- Set order of calculation - if ( td->hasParam("Order") ) { // Local order - hf_errlog(17094510,"W: Ignoring key 'Order' in .dat file. Only global parameter 'Order' is used."); - } - string order = td->getParamS("Order"); // Global order - bool success=true; - // CIJET default is 'NLO' - if (order=="NNLO" ) fnlo->setorder(1); // swith NLO ON - else if (order=="NLO" ) fnlo->setorder(1); // swith NLO ON - else if (order=="LO" ) fnlo->setorder(0); // switch NLO OFF - else hf_errlog(17094502,"E: CIJET. Unrecognized order: "+order); - if (!success) hf_errlog(17094503,"W: CIJET. Requested order "+order+" cannot be set."); - - // --- Set CI normalization - double nm=1.0; - if ( td->hasParam("CInorm") ) { // Local order - nm = *td->getParamD("CInorm"); - } - fnlo->setnorm(nm); - - // --- Set scale factors - double cmur=1, cmuf=1; - if ( td->hasParam("ScaleFacMuR") ) { - hf_errlog(17094504,"I: Setting CIJET scale factor mu_R: "+td->getParamS("ScaleFacMuR")); - cmur = *td->getParamD("ScaleFacMuR"); - } - //else { - // cmur=GetParam("ScaleFacMuR"); - //} - if ( td->hasParam("ScaleFacMuF") ) { - hf_errlog(17094505,"I: Setting CIJET scale factor mu_F: "+td->getParamS("ScaleFacMuF")); - cmuf = *td->getParamD("ScaleFacMuF"); - } - //else { - // cmuf=GetParam("ScaleFacMuF"); - //} - fnlo->setscales(cmuf,cmur); - - // ---Read of the table - fnlo->setinit(); - - } + + int ID = td->id; // ID=dataSetID -> updated to termdata td + if ( td->hasParam("Filename") ) { + // --- Read file and instantiate CIJET + const std::string filename = td->getParamS("Filename"); + hf_errlog(17086501,"I: Reading CIJET file "+filename); + ffnlo.insert(std::make_pair(ID,vector<CIJETReaction*>{new CIJETReaction(filename,this)} )); + } + else if ( td->hasParam("Filenames") ) { + const std::string filenames = td->getParamS("Filenames"); + std::stringstream ss(filenames); + std::string filename; + while (std::getline(ss, filename, ',')) { + hf_errlog(17086501,"I: Reading CIJET file "+filename); + ffnlo[ID].push_back(new CIJETReaction(filename,this)); + } + } + else { + hf_errlog(17086503,"F:No CIJET file specified. Please provide 'Filename' or 'Filenames' to reaction CIJET."); + return; + } + + for ( CIJETReaction* fnlo : ffnlo[ID] ) { + // --- Set order of calculation + if ( td->hasParam("Order") ) { // Local order + hf_errlog(17094510,"W: Ignoring key 'Order' in .dat file. Only global parameter 'Order' is used."); + } + string order = td->getParamS("Order"); // Global order + bool success=true; + // CIJET default is 'NLO' + if (order=="NNLO") fnlo->setorder(1); // switch NLO ON + else if (order=="NLO" ) fnlo->setorder(1); // switch NLO ON + else if (order=="LO" ) fnlo->setorder(0); // switch NLO OFF + else hf_errlog(17094502,"E: CIJET. Unrecognized order: "+order); + if (!success) hf_errlog(17094503,"W: CIJET. Requested order "+order+" cannot be set."); + + // --- Set CI normalization + double nm=1.0; + if ( td->hasParam("CInorm") ) nm = *td->getParamD("CInorm"); // Local order + fnlo->setnorm(nm); + + // --- Set scale factors + double cmur=1, cmuf=1; + if ( td->hasParam("ScaleFacMuR") ) { + hf_errlog(17094504,"I: Setting CIJET scale factor mu_R: "+td->getParamS("ScaleFacMuR")); + cmur = *td->getParamD("ScaleFacMuR"); + } + if ( td->hasParam("ScaleFacMuF") ) { + hf_errlog(17094505,"I: Setting CIJET scale factor mu_F: "+td->getParamS("ScaleFacMuF")); + cmuf = *td->getParamD("ScaleFacMuF"); + } + fnlo->setscales(cmuf,cmur); + + // ---Read the table + fnlo->setinit(); + + } } //______________________________________________________________________________ // Main function to compute results at an iteration void ReactionCIJET::compute(TermData* td, valarray<double> &val, map<string, valarray<double> > &err) { - td->actualizeWrappers(); - // Get relevant parameters: - int cicase = td->getParamI("CIcase"); - double lamotev = *td->getParamD("LamoTeV"); - double cp1 = *td->getParamD("CI1"); - double cp3 = *td->getParamD("CI3"); - double cp5 = *td->getParamD("CI5"); - // Adjust according to type of contact interactions - // 1-3 constrained fit - if ( cicase==1 ) { //pure left or right handed - cp3=0.0; //using cp1 as free parameter - cp5=0.0; - } - else if ( cicase==2 ) { //vector case - cp3=2*cp1; //using cp1 as free parameter - cp5=cp1; - } - else if ( cicase==3 ) { //axial-vector case - cp3=-2*cp1; //using cp1 as free parameter - cp5=cp1; - } - else { //default free-N parameters - } - - vector<double> ci; - ci.push_back(cp1); - ci.push_back(cp3); - ci.push_back(cp5); - - vector<double> valfnlo; - valfnlo.reserve(val.size()); - int ID = td->id; - for ( CIJETReaction* fnlo : ffnlo[ID] ) { - fnlo->setci(lamotev,ci); - std::vector<double> cs = fnlo->calcxsec(); - for ( std::size_t i=0; i<cs.size(); i++) valfnlo.push_back(cs[i]); - } - - if ( val.size()!=valfnlo.size() ) - hf_errlog(17116501,"F: Size of CIJET cross section array does not match data."); - for ( std::size_t i=0; i<val.size(); i++) val[i]=(i<valfnlo.size())?valfnlo[i]:0.; + td->actualizeWrappers(); + // Get relevant parameters: + int cicase = td->getParamI("CIcase"); + double lamotev = *td->getParamD("LamoTeV"); + double cp1 = *td->getParamD("CI1"); + double cp3 = *td->getParamD("CI3"); + double cp5 = *td->getParamD("CI5"); + // Adjust according to type of contact interactions + // 1-3 constrained fit + if ( cicase==1 ) { //pure left or right handed + cp3=0.0; //using cp1 as free parameter + cp5=0.0; + } + else if ( cicase==2 ) { //vector case + cp3=2*cp1; //using cp1 as free parameter + cp5=cp1; + } + else if ( cicase==3 ) { //axial-vector case + cp3=-2*cp1; //using cp1 as free parameter + cp5=cp1; + } + else { //default free-N parameters + } + + vector<double> ci; + ci.push_back(cp1); + ci.push_back(cp3); + ci.push_back(cp5); + + vector<double> valfnlo; + valfnlo.reserve(val.size()); + int ID = td->id; + for ( CIJETReaction* fnlo : ffnlo[ID] ) { + fnlo->setci(lamotev,ci); + std::vector<double> cs = fnlo->calcxsec(); + for ( std::size_t i=0; i<cs.size(); i++) valfnlo.push_back(cs[i]); + } + + if ( val.size()!=valfnlo.size() ) + hf_errlog(17116501,"F: Size of CIJET cross section array does not match data."); + for ( std::size_t i=0; i<val.size(); i++) val[i]=(i<valfnlo.size())?valfnlo[i]:0.; } - - - diff --git a/reactions/CIJET/src/cijet.f b/reactions/CIJET/src/cijet.f index cbe8f50ae..ba20340a3 100644 --- a/reactions/CIJET/src/cijet.f +++ b/reactions/CIJET/src/cijet.f @@ -6,12 +6,11 @@ c performing grid reading and initialization ! new line implicit none ! maximum number of data sets and points per set - integer, parameter :: mset=1, mpint=200 + integer, parameter :: mset=4, mpint=200 integer, parameter :: xnd=40, qnd=15, nbk=42 ! input scale choices, CI shape choice, only color-singlet real(8) :: mufip, murip, fut - integer :: cshape ! 1 for LL, 2 for RR, 3 for vector, 4 for av real(8) :: cpl(3)=0.d0 ! color singlet couplings, LL, LR, RR real(8) :: ash(10) ! individual coupling combinations ! input QCD order @@ -75,15 +74,7 @@ c performing grid reading and initialization ! end of initialization of commons -! define the couplings apart from 1/Lambda^2 -! select case (cshape) -! case(3) -! cpl(1)=1.d0;cpl(2)= 2.d0;cpl(3)=1.d0 -! case(4) -! cpl(1)=1.d0;cpl(2)=-2.d0;cpl(3)=1.d0 -! case default -! cpl(1)=1.d0;cpl(2)=0.d0;cpl(3)=0.d0 -! end select +! init couplings apart from 1/Lambda^2 ash(1) =1.d0 !cpl(1)+cpl(3) ash(2) =1.d0 !cpl(1)^2+cpl(3)^2 ash(3) =1.d0 !cpl(2)^2 @@ -314,7 +305,7 @@ c perform grid convolution and calculation of xsecs ! new line implicit none ! maximum number of data sets and points per set - integer, parameter :: mset=1, mpint=200 + integer, parameter :: mset=4, mpint=200 integer, parameter :: xnd=40, qnd=15, nbk=42 ! input 1/Lambda^2 [in TeV^-2] real(8) :: invlamsq, res(mpint), pres(mpint) @@ -514,9 +505,7 @@ c end of PDF calculation end subroutine cijetxsec -!!!!!!must rely on external Fortran routines for alphas -!!!!!!and PDFs; can not really wrapped them into reaction -!!!!!!class due to problem of passing functions within class +!!!!!!Routines for alphas and PDFs c link to PDF routines subroutine cievolPDF(x, q, res) @@ -528,7 +517,7 @@ c link to PDF routines res=0.d0 ! e.g., - call appl_fnpdf(x, Q, xf) + call pdf_xfxq_wrapper(x, Q, xf) ! test with LHAPDF call evolvePDF(x, Q, xf) res(-6:6)=xf @@ -542,12 +531,12 @@ c link to PDF routines ! new line implicit none real(8) :: q - real(8), external :: appl_fnalphas + real(8), external :: alphas_wrapper ! test with LHAPDF real(8), external :: alphasPDF cialphas=0.118d0 ! e.g., - cialphas=appl_fnalphas(Q) + cialphas=alphas_wrapper(Q) ! test with LHAPDF cialphas=alphasPDF(Q) return diff --git a/reactions/CIJET/yaml/parameters_fixed.yaml b/reactions/CIJET/yaml/parameters_fixed.yaml deleted file mode 100644 index b4e990246..000000000 --- a/reactions/CIJET/yaml/parameters_fixed.yaml +++ /dev/null @@ -1,6 +0,0 @@ -CIJET : - LamoTeV : 9.0 - CIcase : 1 -CI1 : 0.00 -CI3 : 0.00 -CI5 : 0.00 diff --git a/reactions/CIJET/yaml/parameters_orig.yaml b/reactions/CIJET/yaml/parameters_orig.yaml deleted file mode 100644 index 0b297fe3a..000000000 --- a/reactions/CIJET/yaml/parameters_orig.yaml +++ /dev/null @@ -1,8 +0,0 @@ -CIJET : - LamoTeV : 10.0 - CIcase : 1 -CI1 : - value: 1.00 - step: 0.01 -CI3 : 1.00 -CI5 : 1.00 -- GitLab From 50c7ff63fa74391d66abf554bd89eb9804264cdd Mon Sep 17 00:00:00 2001 From: Toni Makela <toni.makela@cern.ch> Date: Sun, 21 Nov 2021 14:47:57 +0100 Subject: [PATCH 3/3] Updated CIJET reaction according to requests for merge: directory structure, added check for cset <= mset, things common to cijetinit and cijetxsec etc are now stored in cijetInclude.h --- reactions/CIJET/CIJET.yaml | 6 ++ reactions/CIJET/{src => }/CIJETReader.cc | 8 +- reactions/CIJET/{include => }/CIJETReader.h | 1 + reactions/CIJET/CMakeLists.txt | 10 +-- reactions/CIJET/{src => }/ReactionCIJET.cc | 0 reactions/CIJET/{include => }/ReactionCIJET.h | 1 - reactions/CIJET/{src => }/cijet.f | 73 +++++-------------- reactions/CIJET/cijetInclude.h | 27 +++++++ reactions/CIJET/yaml/parameters.yaml | 6 -- 9 files changed, 63 insertions(+), 69 deletions(-) create mode 100644 reactions/CIJET/CIJET.yaml rename reactions/CIJET/{src => }/CIJETReader.cc (83%) rename reactions/CIJET/{include => }/CIJETReader.h (98%) rename reactions/CIJET/{src => }/ReactionCIJET.cc (100%) rename reactions/CIJET/{include => }/ReactionCIJET.h (97%) rename reactions/CIJET/{src => }/cijet.f (88%) create mode 100644 reactions/CIJET/cijetInclude.h delete mode 100644 reactions/CIJET/yaml/parameters.yaml diff --git a/reactions/CIJET/CIJET.yaml b/reactions/CIJET/CIJET.yaml new file mode 100644 index 000000000..186cfffc4 --- /dev/null +++ b/reactions/CIJET/CIJET.yaml @@ -0,0 +1,6 @@ +#CIJET : + LamoTeV : 10.0 + CIcase : 1 + CI1 : -0.50 + CI3 : 0.00 + CI5 : 0.00 diff --git a/reactions/CIJET/src/CIJETReader.cc b/reactions/CIJET/CIJETReader.cc similarity index 83% rename from reactions/CIJET/src/CIJETReader.cc rename to reactions/CIJET/CIJETReader.cc index 57eb691fa..939ba610d 100644 --- a/reactions/CIJET/src/CIJETReader.cc +++ b/reactions/CIJET/CIJETReader.cc @@ -8,7 +8,7 @@ using namespace std; // Interface to xFitter FORTRAN routines extern "C" { - void cijetinit_(const char* fname, double* mufip, double* murip, int* od, double* norm); + void cijetinit_(const char* fname, double* mufip, double* murip, int* od, double* norm, int* estat); void cijetxsec_(const char* fname, double* lam, double* cl, int* npt, double* pres); } @@ -19,7 +19,11 @@ void CIJETReader::setinit() //note the character length in Fortran is hardwired to 100 char ff[101]; strncpy(ff, filepds.c_str(), 100); - cijetinit_(ff, &muf, &mur, &Oqcd, &fut); + int estat = 0; //Error status for checking amount of CI grids + cijetinit_(ff, &muf, &mur, &Oqcd, &fut, &estat); + if (estat == 21111801) { + hf_errlog(21111801,"F: cset in cijet.f exceeds mset in cijetInclude.h. Increase mset and recompile."); + } } vector<double> CIJETReader::calcxsec(void) diff --git a/reactions/CIJET/include/CIJETReader.h b/reactions/CIJET/CIJETReader.h similarity index 98% rename from reactions/CIJET/include/CIJETReader.h rename to reactions/CIJET/CIJETReader.h index b71fa717f..a210f5489 100644 --- a/reactions/CIJET/include/CIJETReader.h +++ b/reactions/CIJET/CIJETReader.h @@ -12,6 +12,7 @@ #include <fstream> #include <math.h> #include <cstdlib> +#include "hf_errlog.h" #include "ReactionTheory.h" #include "TermData.h" diff --git a/reactions/CIJET/CMakeLists.txt b/reactions/CIJET/CMakeLists.txt index 5df9e3be2..6d277501c 100644 --- a/reactions/CIJET/CMakeLists.txt +++ b/reactions/CIJET/CMakeLists.txt @@ -1,11 +1,9 @@ set(TARGET reactionCIJET) add_library(${TARGET} MODULE - src/ReactionCIJET.cc - src/CIJETReader.cc - src/cijet.f + ReactionCIJET.cc + CIJETReader.cc + cijet.f ) target_link_libraries(${TARGET} PRIVATE xfitter) -target_include_directories(${TARGET} PUBLIC include) install(TARGETS ${TARGET} DESTINATION ${DEST_MODULE}) -install(FILES yaml/parameters.yaml DESTINATION ${DEST_YAML}/reactions RENAME ReactionCIJET.yaml) - +install(FILES CIJET.yaml DESTINATION ${DEST_YAML}/reactions) diff --git a/reactions/CIJET/src/ReactionCIJET.cc b/reactions/CIJET/ReactionCIJET.cc similarity index 100% rename from reactions/CIJET/src/ReactionCIJET.cc rename to reactions/CIJET/ReactionCIJET.cc diff --git a/reactions/CIJET/include/ReactionCIJET.h b/reactions/CIJET/ReactionCIJET.h similarity index 97% rename from reactions/CIJET/include/ReactionCIJET.h rename to reactions/CIJET/ReactionCIJET.h index 80c762541..6168e7a5a 100644 --- a/reactions/CIJET/include/ReactionCIJET.h +++ b/reactions/CIJET/ReactionCIJET.h @@ -8,7 +8,6 @@ //#include <memory> #include <map> #include <vector> -#include "hf_errlog.h" /** @class' ReactionCIJET diff --git a/reactions/CIJET/src/cijet.f b/reactions/CIJET/cijet.f similarity index 88% rename from reactions/CIJET/src/cijet.f rename to reactions/CIJET/cijet.f index ba20340a3..f524f8da3 100644 --- a/reactions/CIJET/src/cijet.f +++ b/reactions/CIJET/cijet.f @@ -1,44 +1,27 @@ -c CIJET1.1 interface made for Xfitter, Jun Gao, 2018.09.10 +c CIJET1.1 interface made for xFitter, Jun Gao, 2018.09.10 +c Updated for 2021 xFitter release, Toni Makela, 2021.11.19 c performing grid reading and initialization - subroutine cijetinit(fname, mufip, murip, oqcd, fut) + subroutine cijetinit(fname, mufip, murip, oqcd, fut, estat) ! new line implicit none + ! maximum number of data sets and points per set - integer, parameter :: mset=4, mpint=200 - integer, parameter :: xnd=40, qnd=15, nbk=42 + include 'cijetInclude.h' -! input scale choices, CI shape choice, only color-singlet +! input scale choices, CI shape choice and couplings, only color-singlet real(8) :: mufip, murip, fut real(8) :: cpl(3)=0.d0 ! color singlet couplings, LL, LR, RR real(8) :: ash(10) ! individual coupling combinations + ! input QCD order integer :: oqcd ! 0 for LO, 1 for NLO -! grid name - character*100 :: fname integer, save :: befirst=1 -! storing total number of sets - integer :: cset - -! storing reference Lambda scale, no of points, fname per set - real(8) :: Nscale(mset) ! hardwired currently - integer :: pset(mset) - character*100 :: ffnm(mset) -! storing pp or ppbar, and mur/muf - integer :: sig(mset) - real(8) :: srof(mset) - -! storing reference x-Q-mu scale per point - real(8) :: xgd(mset,mpint,xnd) - real(8) :: qgd(mset,mpint,qnd) - real(8) :: smu(mset,mpint) - integer :: nxgd(mset,mpint) - integer :: nqgd(mset,mpint) - -! storing grid weight per point - real(8) :: ciwgt(10,xnd,xnd*qnd,5,mset,mpint) + +! store non-zero error status identifying number if need be + integer :: estat ! beta function needed for scale evaluation real(8), parameter :: beta0=(11*3.d0-2*5.d0)/12.d0/3.1415926d0 @@ -120,6 +103,12 @@ c performing grid reading and initialization ! update record cset=cset+1 + if((cset.gt.mset)) then + estat=21111801 + print *, "F: cset in xfitter.f exceeds mset in xfitterInclude.h." + print *, " Increase mset and recompile.",fname; + return + endif pset(cset)=(nt-1)/(3*nbk) ffnm(cset)=fname sig(cset)=1 @@ -304,9 +293,10 @@ c perform grid convolution and calculation of xsecs ! new line implicit none + ! maximum number of data sets and points per set - integer, parameter :: mset=4, mpint=200 - integer, parameter :: xnd=40, qnd=15, nbk=42 + include 'cijetInclude.h' + ! input 1/Lambda^2 [in TeV^-2] real(8) :: invlamsq, res(mpint), pres(mpint) ! sign (+ for destructive, -for con.) @@ -314,29 +304,6 @@ c perform grid convolution and calculation of xsecs real(8) :: ash(10) ! individual coupling combinations ! number of data points and additional multi-factor integer :: npt -! grid name - character*100 :: fname - -! storing total number of sets - integer :: cset - -! storing reference Lambda scale, no of points, fname per set - real(8) :: Nscale(mset) ! hardwired currently - integer :: pset(mset) - character*100 :: ffnm(mset) -! storing pp or ppbar, and mur/muf - integer :: sig(mset) - real(8) :: srof(mset) - -! storing reference x-Q-mu scale per point - real(8) :: xgd(mset,mpint,xnd) - real(8) :: qgd(mset,mpint,qnd) - real(8) :: smu(mset,mpint) - integer :: nxgd(mset,mpint) - integer :: nqgd(mset,mpint) - -! storing grid weight per point - real(8) :: ciwgt(10,xnd,xnd*qnd,5,mset,mpint) ! workspace real(8) :: XPDF(-6:6, xnd, qnd), xi(xnd), qi(qnd) @@ -360,8 +327,6 @@ c perform grid convolution and calculation of xsecs print *, "CI grid not loaded!",fname; stop endif -!test print *, "gaojun1 !", invlamsq, cpl, cset, pid, pset(pid) - ! find the correct couplings ash(1) =cpl(1)+cpl(3) ash(2) =cpl(1)**2+cpl(3)**2 diff --git a/reactions/CIJET/cijetInclude.h b/reactions/CIJET/cijetInclude.h new file mode 100644 index 000000000..838aea719 --- /dev/null +++ b/reactions/CIJET/cijetInclude.h @@ -0,0 +1,27 @@ +! maximum number of data sets and points per set + integer, parameter :: mset=4, mpint=200 + integer, parameter :: xnd=40, qnd=15, nbk=42 + +! grid name + character*100 :: fname + +! storing total number of sets + integer :: cset + +! storing reference Lambda scale, no of points, fname per set + real(8) :: Nscale(mset) ! hardwired currently + integer :: pset(mset) + character*100 :: ffnm(mset) +! storing pp or ppbar, and mur/muf + integer :: sig(mset) + real(8) :: srof(mset) + +! storing reference x-Q-mu scale per point + real(8) :: xgd(mset,mpint,xnd) + real(8) :: qgd(mset,mpint,qnd) + real(8) :: smu(mset,mpint) + integer :: nxgd(mset,mpint) + integer :: nqgd(mset,mpint) + +! storing grid weight per point + real(8) :: ciwgt(10,xnd,xnd*qnd,5,mset,mpint) diff --git a/reactions/CIJET/yaml/parameters.yaml b/reactions/CIJET/yaml/parameters.yaml deleted file mode 100644 index c707be4be..000000000 --- a/reactions/CIJET/yaml/parameters.yaml +++ /dev/null @@ -1,6 +0,0 @@ -CIJET : - LamoTeV : 10.0 - CIcase : 1 -CI1 : -0.50 -CI3 : 0.00 -CI5 : 0.00 -- GitLab