diff --git a/CMakeLists.txt b/CMakeLists.txt index 0266038e62ae499759d8dcda2479790291350eef..feab2aa94bbb784158becee7514c6e19b964ed11 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -153,6 +153,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/CIJET.yaml b/reactions/CIJET/CIJET.yaml new file mode 100644 index 0000000000000000000000000000000000000000..186cfffc44b7bbf3f40efce1a6b7408c019414fb --- /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/CIJETReader.cc b/reactions/CIJET/CIJETReader.cc new file mode 100644 index 0000000000000000000000000000000000000000..939ba610db5854535bde0ffe96ebd2b4057c4921 --- /dev/null +++ b/reactions/CIJET/CIJETReader.cc @@ -0,0 +1,75 @@ +#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, int* estat); + void cijetxsec_(const char* fname, double* lam, double* cl, int* npt, double* pres); +} + +//------------------------------------------------------------------------------------ + +void CIJETReader::setinit() +{ + //note the character length in Fortran is hardwired to 100 + char ff[101]; + strncpy(ff, filepds.c_str(), 100); + 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) +{ + 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++) 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 +} diff --git a/reactions/CIJET/CIJETReader.h b/reactions/CIJET/CIJETReader.h new file mode 100644 index 0000000000000000000000000000000000000000..a210f5489c814aab15e56cd71f79d42868330fa9 --- /dev/null +++ b/reactions/CIJET/CIJETReader.h @@ -0,0 +1,69 @@ +#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 "hf_errlog.h" +#include "ReactionTheory.h" +#include "TermData.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 + 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; + +}; + +#endif diff --git a/reactions/CIJET/CMakeLists.txt b/reactions/CIJET/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..6d277501c70a6b9b2054c96cf48c8abb582221d2 --- /dev/null +++ b/reactions/CIJET/CMakeLists.txt @@ -0,0 +1,9 @@ +set(TARGET reactionCIJET) +add_library(${TARGET} MODULE + ReactionCIJET.cc + CIJETReader.cc + cijet.f +) +target_link_libraries(${TARGET} PRIVATE xfitter) +install(TARGETS ${TARGET} DESTINATION ${DEST_MODULE}) +install(FILES CIJET.yaml DESTINATION ${DEST_YAML}/reactions) diff --git a/reactions/CIJET/ReactionCIJET.cc b/reactions/CIJET/ReactionCIJET.cc new file mode 100644 index 0000000000000000000000000000000000000000..8fc0e03986813ae74ff04c9974c14a4729a49645 --- /dev/null +++ b/reactions/CIJET/ReactionCIJET.cc @@ -0,0 +1,124 @@ +// 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); // 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.; +} diff --git a/reactions/CIJET/ReactionCIJET.h b/reactions/CIJET/ReactionCIJET.h new file mode 100644 index 0000000000000000000000000000000000000000..6168e7a5a0a56076bc05449a8723b678b82266ea --- /dev/null +++ b/reactions/CIJET/ReactionCIJET.h @@ -0,0 +1,46 @@ +#ifndef xFitter_ReactionCIJET +#define xFitter_ReactionCIJET + +#pragma once + +#include "ReactionTheory.h" +#include "CIJETReader.h" +//#include <memory> +#include <map> +#include <vector> + +/** + @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/cijet.f b/reactions/CIJET/cijet.f new file mode 100644 index 0000000000000000000000000000000000000000..f524f8da398e75d25519d1f394add7338a49bb99 --- /dev/null +++ b/reactions/CIJET/cijet.f @@ -0,0 +1,508 @@ +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, estat) + +! new line + implicit none + +! maximum number of data sets and points per set + include 'cijetInclude.h' + +! 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 + + integer, save :: befirst=1 + +! 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 + 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 + +! 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 + 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 + 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 + 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 + include 'cijetInclude.h' + +! 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 + +! 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 + +! 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 + + +!!!!!!Routines for alphas and PDFs + +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 pdf_xfxq_wrapper(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 :: alphas_wrapper +! test with LHAPDF real(8), external :: alphasPDF + + cialphas=0.118d0 +! e.g., + cialphas=alphas_wrapper(Q) +! test with LHAPDF cialphas=alphasPDF(Q) + return + + end function cialphas diff --git a/reactions/CIJET/cijetInclude.h b/reactions/CIJET/cijetInclude.h new file mode 100644 index 0000000000000000000000000000000000000000..838aea719240af5ed26fccd0bb621d9ebeb977b3 --- /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)