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)