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