From d2c0004f836ac7be035188280d688cacdde13c01 Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Wed, 25 Apr 2018 14:08:09 +0200
Subject: [PATCH 01/17] charm CC implemented in ZMVFNS

---
 .../BaseDISCC/include/ReactionBaseDISCC.h     |   6 +-
 reactions/BaseDISCC/src/ReactionBaseDISCC.cc  | 115 +++++++++++++++---
 2 files changed, 100 insertions(+), 21 deletions(-)

diff --git a/reactions/BaseDISCC/include/ReactionBaseDISCC.h b/reactions/BaseDISCC/include/ReactionBaseDISCC.h
index 4bb3e75d8..fe989d759 100644
--- a/reactions/BaseDISCC/include/ReactionBaseDISCC.h
+++ b/reactions/BaseDISCC/include/ReactionBaseDISCC.h
@@ -34,6 +34,8 @@ class ReactionBaseDISCC : public ReactionTheory
     
     virtual int compute(int dataSetID, valarray<double> &val, map<string, valarray<double> > &err);
   protected:
+    enum class dataFlav { incl, c} ;      //!< Define final state.
+
     virtual int parseOptions(){ return 0;};
 
     virtual void F2 BASE_PARS;
@@ -45,12 +47,14 @@ class ReactionBaseDISCC : public ReactionTheory
     map <int, double>  _polarisation ;      //!< longitudinal polarisation
     map <int, double>  _charge;             //!< lepton beam charge
     map <int, int>    _isReduced;          //!< reduced cross section
+    map <int, dataFlav> _dataFlav;          //!< flavour (incl, c, b)
   protected:
     const int GetNpoint(int dataSetID) {return _npoints[dataSetID];}
     const double GetPolarisation (int dataSetID) {return _polarisation[dataSetID];}
     const double GetCharge(int dataSetID) {return _charge[dataSetID]; }
     const int IsReduced(int dataSetID){ return _isReduced[dataSetID] > 0; }
-    
+    const dataFlav GetDataFlav(int dataSetID) {return _dataFlav[dataSetID]; }
+
     // Another decomposition:
     virtual void GetF2u( int dataSetID, valarray<double>& f2u);
     virtual void GetFLu( int dataSetID, valarray<double>& flu);
diff --git a/reactions/BaseDISCC/src/ReactionBaseDISCC.cc b/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
index 22861d6de..1b6fab5f7 100644
--- a/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
+++ b/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
@@ -10,12 +10,21 @@
 #include <iostream>
 
 // Helpers for QCDNUM (CC): 
+
+//! full
 const double  CCEP2F[] = {0.,0.,1.,0.,1.,0.,0.,1.,0.,1.,0.,0.,0.} ; 
 const double  CCEM2F[] = {0.,0.,0.,1.,0.,1.,0.,0.,1.,0.,1.,0.,0.} ;
 
 const double  CCEP3F[] = {0.,0.,-1.,0.,-1.,0.,0.,1.,0.,1.,0.,0.,0.};
 const double  CCEM3F[] = {0.,0. ,0.,-1.,0.,-1.,0.,0.,1.,0.,1.,0.,0.};
 
+//! c
+const double  CCEP2Fc[] = {0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.} ;
+const double  CCEM2Fc[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.} ;
+
+const double  CCEP3Fc[] = {0.,0.,-1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
+const double  CCEM3Fc[] = {0.,0. ,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.};
+
 // define QCDNUM function:
 extern "C" {
   void zmstfun_(const int& id, const double& key, double& x, double& q2, double& sf, const int& np, const int &flag);
@@ -114,12 +123,21 @@ void  ReactionBaseDISCC::setDatasetParameters( int dataSetID, map<string,string>
 
   // type: sigred, signonred (no F2, FL implemented so far, thus type is defined by bool _isReduced)
   // HERA data files provide 'signonred' CC cross sections
+  // Inclusive "non-reduced" cross section by default.
+  _dataFlav[dataSetID] = dataFlav::incl;
+  string msg = "I: Calculating DIS CC reduced cross section";
   map<string,string>::iterator it = pars.find("type");
   if ( it != pars.end() ) {
-    if(it->second == "sigred")
+    if(it->second == "sigred") 
+    {
       _isReduced[dataSetID] = 1;
+      msg = "I: Calculating DIS CC reduced cross section";
+    }
     else if(it->second == "signonred")
+    {
       _isReduced[dataSetID] = 0;
+      msg = "I: Calculating DIS CC non-reduced cross section";
+    }
     else
     {
       char buffer[256];
@@ -129,6 +147,37 @@ void  ReactionBaseDISCC::setDatasetParameters( int dataSetID, map<string,string>
     }
   }
 
+  // flav: incl, c, b
+  it = pars.find("flav");
+  if ( it != pars.end() ) {
+    if(it->second == "incl")
+    {
+      _dataFlav[dataSetID] = dataFlav::incl;
+      msg += " inclusive";
+    }
+    else if(it->second == "c")
+    {
+      _dataFlav[dataSetID] = dataFlav::c;
+      msg += " charm";
+    }
+    // no beauty
+    else if(it->second == "b")
+    {
+      char buffer[256];
+      sprintf(buffer, "F: predictions for beauty in CC are not available (dataset id = %d)", dataSetID);
+      string str = buffer;
+      hf_errlog_(18042501, str.c_str(), str.length());
+    }
+    else
+    {
+      char buffer[256];
+      sprintf(buffer, "F: dataset with id = %d has unknown flav = %s", dataSetID, it->second.c_str());
+      string str = buffer;
+      hf_errlog_(18042502, str.c_str(), str.length());
+    }
+  }
+  hf_errlog_(17041001, msg.c_str(), msg.size());
+
   // e charge: double
   it = pars.find("echarge");
   if ( it != pars.end() )
@@ -200,11 +249,17 @@ void ReactionBaseDISCC::GetF2u(int dataSetID, valarray<double>& f2u)
     
   // Call QCDNUM
     const int id = 2; const int flag = 0; int Npnt = GetNpoint(dataSetID);
-    zmstfun_(id,CCEP2F[0], x[0], q2[0], (_f2u[dataSetID])[0], Npnt, flag);
-    //  zmstfun_(id,CCEM2F[0], x[0], q2[0], (_f2d[dataSetID])[0], Npnt, flag);    
+    switch ( GetDataFlav(dataSetID) )
+      {
+      case dataFlav::incl :
+        zmstfun_(id,CCEP2F[0], x[0], q2[0], (_f2u[dataSetID])[0], Npnt, flag);
+        break;
+      case dataFlav::c :
+        zmstfun_(id,CCEP2Fc[0], x[0], q2[0], (_f2u[dataSetID])[0], Npnt, flag);
+        break ;
+      }
   }
   f2u = _f2u[dataSetID];
-  //  f2d = _f2d[dataSetID];
 }
 
 void ReactionBaseDISCC::GetFLu(int dataSetID, valarray<double>& flu)
@@ -217,11 +272,17 @@ void ReactionBaseDISCC::GetFLu(int dataSetID, valarray<double>& flu)
     
     // Call QCDNUM
     const int id = 1; const int flag = 0; int Npnt = GetNpoint(dataSetID);
-    zmstfun_(id,CCEP2F[0], x[0], q2[0], (_flu[dataSetID])[0], Npnt, flag);
-    //    zmstfun_(id,CCEM2F[0], x[0], q2[0], (_fld[dataSetID])[0], Npnt, flag);    
+    switch ( GetDataFlav(dataSetID) )
+      {
+      case dataFlav::incl :
+        zmstfun_(id,CCEP2F[0], x[0], q2[0], (_flu[dataSetID])[0], Npnt, flag);
+        break;
+      case dataFlav::c :
+        zmstfun_(id,CCEP2Fc[0], x[0], q2[0], (_flu[dataSetID])[0], Npnt, flag);
+        break ;
+      }
   }
   flu = _flu[dataSetID];
-  // fld = _fld[dataSetID];
 
 }
 
@@ -235,12 +296,17 @@ void ReactionBaseDISCC::GetxF3u( int dataSetID, valarray<double>& xf3u )
     
     // Call QCDNUM
     const int id = 3; const int flag = 0; int Npnt = GetNpoint(dataSetID);
-
-    zmstfun_(id,CCEP3F[0], x[0], q2[0], (_xf3u[dataSetID])[0], Npnt, flag);
-    //    zmstfun_(id,CCEM3F[0], x[0], q2[0], (_xf3d[dataSetID])[0], Npnt, flag);    
+    switch ( GetDataFlav(dataSetID) )
+      {
+      case dataFlav::incl :
+        zmstfun_(id,CCEP3F[0], x[0], q2[0], (_xf3u[dataSetID])[0], Npnt, flag);
+        break;
+      case dataFlav::c :
+        zmstfun_(id,CCEP3Fc[0], x[0], q2[0], (_xf3u[dataSetID])[0], Npnt, flag);
+        break;
+      }
   }
   xf3u = _xf3u[dataSetID];
-  //xf3d = _xf3d[dataSetID];
 }
 
 
@@ -254,10 +320,16 @@ void ReactionBaseDISCC::GetF2d(int dataSetID, valarray<double>& f2d)
     
   // Call QCDNUM
     const int id = 2; const int flag = 0; int Npnt = GetNpoint(dataSetID);
-    // zmstfun_(id,CCEP2F[0], x[0], q2[0], (_f2u[dataSetID])[0], Npnt, flag);
-    zmstfun_(id,CCEM2F[0], x[0], q2[0], (_f2d[dataSetID])[0], Npnt, flag);    
+    switch ( GetDataFlav(dataSetID) )
+      {
+      case dataFlav::incl :
+        zmstfun_(id,CCEM2F[0], x[0], q2[0], (_f2d[dataSetID])[0], Npnt, flag);
+        break;
+      case dataFlav::c :
+        zmstfun_(id,CCEM2Fc[0], x[0], q2[0], (_f2d[dataSetID])[0], Npnt, flag);
+        break ;
+      }
   }
-  //f2u = _f2u[dataSetID];
   f2d = _f2d[dataSetID];
 }
 
@@ -271,10 +343,16 @@ void ReactionBaseDISCC::GetFLd(int dataSetID, valarray<double>& fld)
     
     // Call QCDNUM
     const int id = 1; const int flag = 0; int Npnt = GetNpoint(dataSetID);
-    // zmstfun_(id,CCEP2F[0], x[0], q2[0], (_flu[dataSetID])[0], Npnt, flag);
-    zmstfun_(id,CCEM2F[0], x[0], q2[0], (_fld[dataSetID])[0], Npnt, flag);    
+    switch ( GetDataFlav(dataSetID) )
+      {
+      case dataFlav::incl :
+        zmstfun_(id,CCEM2F[0], x[0], q2[0], (_fld[dataSetID])[0], Npnt, flag);
+        break;
+      case dataFlav::c :
+        zmstfun_(id,CCEM2Fc[0], x[0], q2[0], (_fld[dataSetID])[0], Npnt, flag);
+        break ;
+      }
   }
-  // flu = _flu[dataSetID];
   fld = _fld[dataSetID];
 
 }
@@ -289,11 +367,8 @@ void ReactionBaseDISCC::GetxF3d( int dataSetID, valarray<double>& xf3d )
     
     // Call QCDNUM
     const int id = 3; const int flag = 0; int Npnt = GetNpoint(dataSetID);
-
-    // zmstfun_(id,CCEP3F[0], x[0], q2[0], (_xf3u[dataSetID])[0], Npnt, flag);
     zmstfun_(id,CCEM3F[0], x[0], q2[0], (_xf3d[dataSetID])[0], Npnt, flag);    
   }
-  //xf3u = _xf3u[dataSetID];
   xf3d = _xf3d[dataSetID];
 }
 
-- 
GitLab


From eeb20cf6f6219e0dcbf831d5b74efc08f2aba297 Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Wed, 25 Apr 2018 20:13:50 +0200
Subject: [PATCH 02/17] charm CC in FONLL

---
 .../FONLL_DISCC/src/ReactionFONLL_DISCC.cc    | 20 +++++++++++++++----
 1 file changed, 16 insertions(+), 4 deletions(-)

diff --git a/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc b/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
index c55086271..4cefbaf70 100644
--- a/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
+++ b/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
@@ -178,10 +178,22 @@ void ReactionFONLL_DISCC::initAtIteration()
 	      APFEL::ComputeStructureFunctionsAPFEL(Q,Q);
 	    }
 
-	  // Compute structure functions by interpolation in x.
-	  _f2fonll[dataSetID][i] = APFEL::F2total(x[i]) / 2;
-	  _flfonll[dataSetID][i] = APFEL::FLtotal(x[i]) / 2;
-	  _f3fonll[dataSetID][i] = APFEL::F3total(x[i]) / 2;
+	  // Compute structure functions by interpolation in x for the
+	  // appropriate component (total, charm, or bottom).
+	  switch (GetDataFlav(dataSetID))
+	    {
+	    case dataFlav::incl:
+	      _f2fonll[dataSetID][i] = APFEL::F2total(x[i]) / 2;
+	      _flfonll[dataSetID][i] = APFEL::FLtotal(x[i]) / 2;
+	      _f3fonll[dataSetID][i] = APFEL::F3total(x[i]) / 2;
+	      break;
+	    case dataFlav::c:
+	      _f2fonll[dataSetID][i] = APFEL::F2charm(x[i]) / 2;
+	      _flfonll[dataSetID][i] = APFEL::FLcharm(x[i]) / 2;
+	      _f3fonll[dataSetID][i] = APFEL::F3charm(x[i]) / 2;
+	      break;
+	    }
+	  printf("%.2e  %.2e  %.2e\n", APFEL::F2total(x[i]), APFEL::F2charm(x[i]), APFEL::F2bottom(x[i]));
 
 	  Q2save = q2[i];
 	}
-- 
GitLab


From 50d50bf34c8d7715931485370090258960c70f41 Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Wed, 25 Apr 2018 20:14:33 +0200
Subject: [PATCH 03/17] charm CC in FFABM

---
 .../FFABM_DISCC/src/ReactionFFABM_DISCC.cc      | 17 ++++++-----------
 1 file changed, 6 insertions(+), 11 deletions(-)

diff --git a/reactions/FFABM_DISCC/src/ReactionFFABM_DISCC.cc b/reactions/FFABM_DISCC/src/ReactionFFABM_DISCC.cc
index faf2708a7..334cd8645 100644
--- a/reactions/FFABM_DISCC/src/ReactionFFABM_DISCC.cc
+++ b/reactions/FFABM_DISCC/src/ReactionFFABM_DISCC.cc
@@ -162,24 +162,19 @@ void ReactionFFABM_DISCC::calcF2FL(int dataSetID) {
             }
 
 
-          /*switch ( GetDataType(dataSetID) )
+          switch ( GetDataFlav(dataSetID) )
             {
-            case dataType::sigred :*/
+            case dataFlav::incl :
               _f2abm[dataSetID][i] = f2 + f2c + f2b;
               _flabm[dataSetID][i] = fl + flc + flb;
               _f3abm[dataSetID][i] = x[i] * (f3 + f3c + f3b);
-              /*break;
-            case dataType::f2c :
+              break;
+            case dataFlav::c :
               _f2abm[dataSetID][i] = f2c;
               _flabm[dataSetID][i] = flc;
               _f3abm[dataSetID][i] = x[i] * f3c;
-              break ;
-            case dataType::f2b :
-              _f2abm[dataSetID][i] = f2b;
-              _flabm[dataSetID][i] = flb;
-              _f3abm[dataSetID][i] = x[i] * f3b;
-              break ;
-            }*/
+              break;
+            }
         }
   }
 }
-- 
GitLab


From 50c6ff2ab0627bd2a74ab5f2cb714beaa66226a6 Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Wed, 25 Apr 2018 20:16:39 +0200
Subject: [PATCH 04/17] if y bins not provided, calculate y = Q2 / sx

---
 include/ReactionTheory.h                     |  4 +++-
 reactions/BaseDISCC/src/ReactionBaseDISCC.cc | 17 ++++++++++++++++-
 reactions/BaseDISNC/src/ReactionBaseDISNC.cc | 15 +++++++++++++++
 3 files changed, 34 insertions(+), 2 deletions(-)

diff --git a/include/ReactionTheory.h b/include/ReactionTheory.h
index 8bcb1eeda..3f17cbf2c 100644
--- a/include/ReactionTheory.h
+++ b/include/ReactionTheory.h
@@ -76,7 +76,6 @@ class ReactionTheory
   
   virtual void setBinning(int dataSetID, map<string,valarray<double> > *dsBins){ _dsIDs.push_back(dataSetID); _dsBins[dataSetID] = dsBins; } ;
 
-
   //! Perform optional re-initialization for a given iteration
   virtual void initAtIteration() {};
 
@@ -193,6 +192,9 @@ class ReactionTheory
     }
   };
 
+  // Add one more array of bins which can be calculated using another provided arrays, but suitable to store and use in caluclations (e.g. in DIS y = Q2 / sx)
+  virtual void AddBinning(int dataSetID, std::pair<string,valarray<double>* >* dsBin){ (*_dsBins[dataSetID])[dsBin->first] = *dsBin->second; } ;
+
  protected:
   string _subtype;
   valarray<double> *_val;
diff --git a/reactions/BaseDISCC/src/ReactionBaseDISCC.cc b/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
index 1b6fab5f7..07e07be08 100644
--- a/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
+++ b/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
@@ -109,7 +109,22 @@ void ReactionBaseDISCC::initAtIteration() {
 // 
 void  ReactionBaseDISCC::setDatasetParameters( int dataSetID, map<string,string> pars, map<string,double> parsDataset) 
 {
-  auto *q2p  = GetBinValues(dataSetID,"Q2"), *xp  = GetBinValues(dataSetID,"x"), *yp  = GetBinValues(dataSetID,"y");  
+  auto *q2p  = GetBinValues(dataSetID,"Q2"), *xp  = GetBinValues(dataSetID,"x"), *yp  = GetBinValues(dataSetID,"y");
+
+  // if Q2 and x bins and centre-of-mass energy provided, calculate y = Q2 / (s * x)
+  if(yp == nullptr && q2p != nullptr && xp != nullptr)
+  {
+    map<string,string>::iterator it = pars.find("energy");
+    if ( it != pars.end() )
+    {
+      double s = pow(stof(it->second), 2.0);
+      valarray<double> y = (*q2p) / (s * (*xp));
+      std::pair<string,valarray<double>* > dsBin = std::make_pair("y", &y);
+      AddBinning(dataSetID, &dsBin);
+      yp = GetBinValues(dataSetID, "y");
+    }
+  }
+
   if (q2p == nullptr || xp == nullptr || yp == nullptr ) {
     string msg = "F: Q2, x or Y bins are missing for CC DIS reaction for dataset " + std::to_string(dataSetID);
     hf_errlog_(17100801,msg.c_str(), msg.size());
diff --git a/reactions/BaseDISNC/src/ReactionBaseDISNC.cc b/reactions/BaseDISNC/src/ReactionBaseDISNC.cc
index 701b6be2b..e26927300 100644
--- a/reactions/BaseDISNC/src/ReactionBaseDISNC.cc
+++ b/reactions/BaseDISNC/src/ReactionBaseDISNC.cc
@@ -94,6 +94,21 @@ void ReactionBaseDISNC::initAtIteration() {
 void  ReactionBaseDISNC::setDatasetParameters( int dataSetID, map<string,string> pars, map<string,double> parsDataset) 
 {
   auto *q2p  = GetBinValues(dataSetID,"Q2"), *xp  = GetBinValues(dataSetID,"x"), *yp  = GetBinValues(dataSetID,"y");  
+
+  // if Q2 and x bins and centre-of-mass energy provided, calculate y = Q2 / (s * x)
+  if(yp == nullptr && q2p != nullptr && xp != nullptr)
+  {
+    map<string,string>::iterator it = pars.find("energy");
+    if ( it != pars.end() )
+    {
+      double s = pow(stof(it->second), 2.0);
+      valarray<double> y = (*q2p) / (s * (*xp));
+      std::pair<string,valarray<double>* > dsBin = std::make_pair("y", &y);
+      AddBinning(dataSetID, &dsBin);
+      yp = GetBinValues(dataSetID, "y");
+    }
+  }
+
   if (q2p == nullptr || xp == nullptr || yp == nullptr ) {
     string msg = "F: Q2, x or Y bins are missing for NC DIS reaction for dataset " + std::to_string(dataSetID);
     hf_errlog_(17040801,msg.c_str(), msg.size());
-- 
GitLab


From 2da6dc3bbb80562d077f475fcd2ebf17f0db1ee2 Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Wed, 25 Apr 2018 20:34:41 +0200
Subject: [PATCH 05/17] disabled printout

---
 reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc | 1 -
 1 file changed, 1 deletion(-)

diff --git a/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc b/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
index 4cefbaf70..b2a656ab6 100644
--- a/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
+++ b/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
@@ -193,7 +193,6 @@ void ReactionFONLL_DISCC::initAtIteration()
 	      _f3fonll[dataSetID][i] = APFEL::F3charm(x[i]) / 2;
 	      break;
 	    }
-	  printf("%.2e  %.2e  %.2e\n", APFEL::F2total(x[i]), APFEL::F2charm(x[i]), APFEL::F2bottom(x[i]));
 
 	  Q2save = q2[i];
 	}
-- 
GitLab


From 6b4b2f397d3ea23c148e100c99f7b64e5462993a Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Wed, 9 May 2018 14:07:49 +0200
Subject: [PATCH 06/17] further implemententation, validation, tests

---
 reactions/BaseDISCC/src/ReactionBaseDISCC.cc | 42 ++++++++++++++++----
 src/xfitter_pars.cc                          |  7 ++++
 2 files changed, 41 insertions(+), 8 deletions(-)

diff --git a/reactions/BaseDISCC/src/ReactionBaseDISCC.cc b/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
index 07e07be08..8430d02b9 100644
--- a/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
+++ b/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
@@ -12,18 +12,36 @@
 // Helpers for QCDNUM (CC): 
 
 //! full
-const double  CCEP2F[] = {0.,0.,1.,0.,1.,0.,0.,1.,0.,1.,0.,0.,0.} ; 
-const double  CCEM2F[] = {0.,0.,0.,1.,0.,1.,0.,0.,1.,0.,1.,0.,0.} ;
+const double  CCEP2F[] = {0.,0.,1.,0.,1.,0., 0. ,1.,0.,1.,0.,0.,0.} ;
+const double  CCEM2F[] = {0.,0.,0.,1.,0.,1., 0. ,0.,1.,0.,1.,0.,0.} ;
 
 const double  CCEP3F[] = {0.,0.,-1.,0.,-1.,0.,0.,1.,0.,1.,0.,0.,0.};
 const double  CCEM3F[] = {0.,0. ,0.,-1.,0.,-1.,0.,0.,1.,0.,1.,0.,0.};
 
 //! c
-const double  CCEP2Fc[] = {0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.} ;
-const double  CCEM2Fc[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.} ;
-
-const double  CCEP3Fc[] = {0.,0.,-1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
-const double  CCEM3Fc[] = {0.,0. ,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.};
+// work in progress: according to 1001.2312 section 5,
+// in ZM only the sum of contributions s + c makes sense
+// three different options are below for checks, uncommented one is for s + c
+//
+// only c
+//const double  CCEP2Fc[] = {0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.} ;
+//const double  CCEM2Fc[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.} ;
+// only s
+//const double  CCEP2Fc[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,0.} ;
+//const double  CCEM2Fc[] = {0.,0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.} ;
+// only s,c
+const double  CCEP2Fc[] = {0.,0.,1.,0.,0.,0.,0.,0.,0.,1.,0.,0.,0.} ;
+const double  CCEM2Fc[] = {0.,0.,0.,1.,0.,0.,0.,0.,0.,0.,1.,0.,0.} ;
+
+// only c
+//const double  CCEP3Fc[] = {0.,0.,-1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
+//const double  CCEM3Fc[] = {0.,0. ,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.};
+// only s
+//const double  CCEP3Fc[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,0.};
+//const double  CCEM3Fc[] = {0.,0.,0.,-1.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
+// only s,c
+const double  CCEP3Fc[] = {0.,0.,-1.,0.,0.,0.,0.,0.,0.,1.,0.,0.,0.};
+const double  CCEM3Fc[] = {0.,0.,0.,-1.,0.,0.,0.,0.,0.,0.,1.,0.,0.};
 
 // define QCDNUM function:
 extern "C" {
@@ -382,7 +400,15 @@ void ReactionBaseDISCC::GetxF3d( int dataSetID, valarray<double>& xf3d )
     
     // Call QCDNUM
     const int id = 3; const int flag = 0; int Npnt = GetNpoint(dataSetID);
-    zmstfun_(id,CCEM3F[0], x[0], q2[0], (_xf3d[dataSetID])[0], Npnt, flag);    
+    switch ( GetDataFlav(dataSetID) )
+      {
+      case dataFlav::incl :
+        zmstfun_(id,CCEM3F[0], x[0], q2[0], (_xf3d[dataSetID])[0], Npnt, flag);
+        break;
+      case dataFlav::c :
+        zmstfun_(id,CCEM3Fc[0], x[0], q2[0], (_xf3d[dataSetID])[0], Npnt, flag);
+        break;
+    }
   }
   xf3d = _xf3d[dataSetID];
 }
diff --git a/src/xfitter_pars.cc b/src/xfitter_pars.cc
index 55f10c412..d3eaf499b 100644
--- a/src/xfitter_pars.cc
+++ b/src/xfitter_pars.cc
@@ -216,6 +216,13 @@ namespace XFITTER_PARS {
     // constants
     FortAssignD(Gf,constants_)
     FortAssignD(ConvFac,constants_)
+    // OZ 26.04.18 the two lines above do not do what expected because it is gf and convFac in parameters.yaml, not Gf and ConvFac
+    // as a result CC DIS cross section with old interface are always zero
+    // temporary fix: set these parameters manually
+    if(gParameters.find("gf") != gParameters.end())
+      constants_.Gf = *gParameters["gf"];
+    if(gParameters.find("convFac") != gParameters.end())
+      constants_.ConvFac = *gParameters["convFac"];
 
     //Fermion masses:
     FortAssignD(men,fermion_masses_)  // electron neutrino
-- 
GitLab


From 67352ae9dd9cba13a5120bb5d3a9e9f8eeb4d8b7 Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Sun, 3 Jun 2018 00:31:14 +0200
Subject: [PATCH 07/17] one more fir for parameter reading

---
 include/IntegrateDIS.h |  0
 src/IntegrateDIS.cc    |  0
 src/xfitter_pars.cc    | 10 +++++++---
 3 files changed, 7 insertions(+), 3 deletions(-)
 create mode 100644 include/IntegrateDIS.h
 create mode 100644 src/IntegrateDIS.cc

diff --git a/include/IntegrateDIS.h b/include/IntegrateDIS.h
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/IntegrateDIS.cc b/src/IntegrateDIS.cc
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/xfitter_pars.cc b/src/xfitter_pars.cc
index d3eaf499b..5ee93133f 100644
--- a/src/xfitter_pars.cc
+++ b/src/xfitter_pars.cc
@@ -216,9 +216,13 @@ namespace XFITTER_PARS {
     // constants
     FortAssignD(Gf,constants_)
     FortAssignD(ConvFac,constants_)
-    // OZ 26.04.18 the two lines above do not do what expected because it is gf and convFac in parameters.yaml, not Gf and ConvFac
-    // as a result CC DIS cross section with old interface are always zero
-    // temporary fix: set these parameters manually
+
+    // OZ 26.04.18 some lines above do not do what expected because it is gf, convFac and alphaem in parameters.yaml, not Gf, ConvFac and Alphaem
+    // as a result CC DIS cross section and integrated NC and CC cross sections are always zero with old interface
+    // temporary fix: set these parameters manually 
+    // (maybe some other parameters are not assigned as well)
+    if(gParameters.find("alphaem") != gParameters.end())
+      ew_couplings_.Alphaem = *gParameters["alphaem"];
     if(gParameters.find("gf") != gParameters.end())
       constants_.Gf = *gParameters["gf"];
     if(gParameters.find("convFac") != gParameters.end())
-- 
GitLab


From 263c345899e6563e4def7c43d2a9d476a885045e Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Sun, 3 Jun 2018 00:32:08 +0200
Subject: [PATCH 08/17] code for integrated DIS

---
 include/IntegrateDIS.h |  41 +++++++++++++
 src/IntegrateDIS.cc    | 129 +++++++++++++++++++++++++++++++++++++++++
 src/Makefile.am        |   2 +
 3 files changed, 172 insertions(+)

diff --git a/include/IntegrateDIS.h b/include/IntegrateDIS.h
index e69de29bb..c7a0b08d7 100644
--- a/include/IntegrateDIS.h
+++ b/include/IntegrateDIS.h
@@ -0,0 +1,41 @@
+#include <valarray>
+
+// Class used by ReactionBaseDISNC and ReactionBaseDISCC for integrating DIS cross sections
+// and providing them over Q2, y, x ranges.
+// Method is based on legacy subroutine GetIntegratedDisXsection from dis_sigma.f
+
+class IntegrateDIS
+{
+  private:
+    // number of Q2, x subbins are hardcoded now
+    // in the future they could be specified optionally in data files and passed here from DIS reactions classes
+    const int _nsplit_x = 25;
+    const int _nsplit_q2 = 25;
+    //const int _nsplit_x = 100;
+    //const int _nsplit_q2 = 100;
+    //const int _nsplit_x = 5;
+    //const int _nsplit_q2 = 5;
+    std::valarray<int> _nSubBins;
+    std::valarray<double> _q2;
+    std::valarray<double> _x;
+    std::valarray<double> _y;
+    std::valarray<double> _deltaq2;
+    std::valarray<double> _deltax;
+
+  public:
+    //IntegrateDIS();
+
+    // initialise integrated cross section for one dataset, return number of subbins
+    int init(const double s,
+             std::valarray<double>* q2minp, std::valarray<double>* q2maxp,
+             std::valarray<double>* yminp, std::valarray<double>* ymaxp,
+             std::valarray<double>* xminp, std::valarray<double>* xmaxp);
+
+    // calculate integrated cross sections by integrating over subbins
+    std::valarray<double> compute(const std::valarray<double>& val);
+
+    // get bin values
+    std::valarray<double>* getBinValuesQ2() { return &_q2; }
+    std::valarray<double>* getBinValuesX() { return &_x; }
+    std::valarray<double>* getBinValuesY() { return &_y; }
+};
diff --git a/src/IntegrateDIS.cc b/src/IntegrateDIS.cc
index e69de29bb..02e1cc64e 100644
--- a/src/IntegrateDIS.cc
+++ b/src/IntegrateDIS.cc
@@ -0,0 +1,129 @@
+#include <IntegrateDIS.h>
+#include <xfitter_cpp_base.h>
+#include <cassert>
+
+int IntegrateDIS::init(const double s,
+                        std::valarray<double>* q2minp, std::valarray<double>* q2maxp,
+                        std::valarray<double>* yminp, std::valarray<double>* ymaxp,
+                        std::valarray<double>* xminp, std::valarray<double>* xmaxp)
+{
+  // prepare number of bins and arrays
+  const int npoints = q2minp->size();
+  _nSubBins.resize(npoints);
+  int nSubBins = _nsplit_x * _nsplit_q2;
+  _q2.resize(npoints * nSubBins);
+  _deltaq2.resize(npoints * nSubBins);
+  _x.resize(npoints * nSubBins);
+  _deltax.resize(npoints * nSubBins);
+  _y.resize(npoints * nSubBins);
+  int currentSubBin = 0;
+
+  for(size_t i = 0; i < q2minp->size(); i++)
+  {
+    const double q2min = (*q2minp)[i];
+    const double q2max = (*q2maxp)[i];
+    const double ymin = (*yminp)[i];
+    const double ymax = (*ymaxp)[i];
+
+    //if(ymin == 0.)
+
+    // just copy previous entry, if binning is the same
+    //bool flagCopyValue = false;
+    if(i > 0 && q2min == (*q2minp)[i - 1] && q2max == (*q2maxp)[i - 1] && ymin == (*yminp)[i - 1] && ymax == (*ymaxp)[i - 1])
+      //flagCopyValue = true;
+      //if(flagCopyValue)
+    {
+      // -1 means that the value from previous bin will be copied in compute()
+      _nSubBins[i] = -1;
+      continue;
+    }
+
+    // determine x range
+    double xmin = q2min / (s * ymax);
+    double xmax = q2max / (s * ymin);
+    if(xminp)
+      xmin = std::max(xmin, (*xminp)[i]);
+    if(xmaxp)
+      xmax = std::min(xmax, (*xmaxp)[i]);
+    if(xmax > 1.0)
+      xmax = 0.999999;
+      //hf_errlog(18060101, "xmaxCalc > 1.0");
+
+    // do integration in log space for Q2 and x
+    int j = 0;
+    double q2_1 = -1.0;
+    double q2_2 = -1.0;
+    double x_1 = -1.0;
+    double x_2 = -1.0;
+    for(int iq2 = 0; iq2 <= _nsplit_q2; iq2++)
+    {
+      q2_1 = q2_2;
+      q2_2 = exp(log(q2min) + (log(q2max) - log(q2min)) / _nsplit_q2 * iq2);
+      if(iq2 > 0)
+      {
+        for(int ix = 0; ix <= _nsplit_x; ix++)
+        {
+          x_1 = x_2;
+          x_2 = exp(log(xmin) + (log(xmax) - log(xmin)) / _nsplit_x * ix);
+          if(ix > 0)
+          {
+            j = j + 1;
+            _deltaq2[currentSubBin] = q2_2 - q2_1;
+            double q2 = exp(log(q2_1) + 0.5 * (log(q2_2) - log(q2_1)));
+            _q2[currentSubBin] = q2;
+            _deltax[currentSubBin] = x_2 - x_1;
+            double x = exp(log(x_1) + 0.5 * (log(x_2) - log(x_1)));
+            _x[currentSubBin] = x;
+            double y = q2 / (s * x);
+            // check that calculated y values agree with limits given in data file
+            // TODO: this is not optimal, better would be to skip such points
+            // however, then arrays of subbins will have different size
+            if(y < ymin || y > ymax)
+              y = 0.0;
+            _y[currentSubBin] = y;
+            currentSubBin++;
+          }
+        }
+      }
+    }
+    assert(nSubBins == j);
+    _nSubBins[i] = nSubBins;
+  }
+  return currentSubBin;
+}
+
+std::valarray<double> IntegrateDIS::compute(const std::valarray<double>& val)
+{
+  std::valarray<double> valIntegrated(_nSubBins.size());
+  unsigned int nSubBins = 0;
+  unsigned int nBins = 0;
+  for(auto& n : _nSubBins)
+  {
+    if(n == -1)
+    {
+      // copy cross section from previous bin
+      valIntegrated[nBins] = valIntegrated[nBins - 1];
+    }
+    else
+    {
+      // integrate
+      double xsec = 0.0;
+      for(size_t i = nSubBins; i < nSubBins + n; i++)
+      {
+        // skip points out of kinematic range
+        if(_y[i] == 0.0)
+          continue;
+        double dxsec = val[i] * _deltaq2[i] * _deltax[i];
+        xsec += dxsec;
+        printf("%f %f: %f += %f [ %f * %f * %f ]\n", _q2[i], _x[i],
+               xsec, dxsec, val[i], _deltaq2[i], _deltax[i]);
+      }
+      valIntegrated[nBins] = xsec;
+      // increase number of processed bins
+      nSubBins += n;
+    }
+    nBins++;
+  }
+  assert(valIntegrated.size() == nBins);
+  return valIntegrated;
+}
diff --git a/src/Makefile.am b/src/Makefile.am
index 310680df1..95466bc6f 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -38,6 +38,8 @@ libxfmain_la_SOURCES += fortran_io.f lhapdferrors.cc
 libxfmain_la_SOURCES += chi2scan.cc
 # ReactionTheory
 libxfmain_la_SOURCES += ReactionTheory.cc
+# IntegrateDIS
+libxfmain_la_SOURCES += IntegrateDIS.cc
 # optional sources
 # these should be moved to indiviual shared libraries in the future
 if ENABLE_APFEL
-- 
GitLab


From 87f5aaf4778ea699a8689ed02e9065c43ade3e43 Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Sun, 3 Jun 2018 00:32:48 +0200
Subject: [PATCH 09/17] made GetBinValues virtual, needed for integrated DIS

---
 include/ReactionTheory.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/ReactionTheory.h b/include/ReactionTheory.h
index 3f17cbf2c..281884365 100644
--- a/include/ReactionTheory.h
+++ b/include/ReactionTheory.h
@@ -175,7 +175,7 @@ class ReactionTheory
 
 
   // Helper function to get bin values for a given data set, bin name. Returns null if not found
-  valarray<double> *GetBinValues(int idDS, const string& binName)
+  virtual valarray<double> *GetBinValues(int idDS, const string& binName)
   { 
     map<string, valarray<double> >* mapBins =  _dsBins[idDS];
     if (mapBins == nullptr ) {
-- 
GitLab


From 206352d7ed5cbac4e171919520c7bda364b859fe Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Sun, 3 Jun 2018 00:33:48 +0200
Subject: [PATCH 10/17] disabled printout

---
 src/IntegrateDIS.cc | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/IntegrateDIS.cc b/src/IntegrateDIS.cc
index 02e1cc64e..8a2a1138e 100644
--- a/src/IntegrateDIS.cc
+++ b/src/IntegrateDIS.cc
@@ -115,8 +115,8 @@ std::valarray<double> IntegrateDIS::compute(const std::valarray<double>& val)
           continue;
         double dxsec = val[i] * _deltaq2[i] * _deltax[i];
         xsec += dxsec;
-        printf("%f %f: %f += %f [ %f * %f * %f ]\n", _q2[i], _x[i],
-               xsec, dxsec, val[i], _deltaq2[i], _deltax[i]);
+        //printf("%f %f: %f += %f [ %f * %f * %f ]\n", _q2[i], _x[i],
+        //       xsec, dxsec, val[i], _deltaq2[i], _deltax[i]);
       }
       valIntegrated[nBins] = xsec;
       // increase number of processed bins
-- 
GitLab


From ac334d02bbcf3eb3a7ec08ddfd354c6ffa14874d Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Sun, 3 Jun 2018 00:36:27 +0200
Subject: [PATCH 11/17] implemented integrated DIS, updated CC charm for ZMVFNS
 to both s+c initial states

---
 .../BaseDISCC/include/ReactionBaseDISCC.h     |  10 +-
 reactions/BaseDISCC/src/Makefile.am           |   2 +
 reactions/BaseDISCC/src/ReactionBaseDISCC.cc  | 136 ++++++++++++----
 .../BaseDISNC/include/ReactionBaseDISNC.h     |  14 +-
 reactions/BaseDISNC/src/ReactionBaseDISNC.cc  | 147 ++++++++++++++----
 5 files changed, 254 insertions(+), 55 deletions(-)

diff --git a/reactions/BaseDISCC/include/ReactionBaseDISCC.h b/reactions/BaseDISCC/include/ReactionBaseDISCC.h
index fe989d759..cfa0fae8b 100644
--- a/reactions/BaseDISCC/include/ReactionBaseDISCC.h
+++ b/reactions/BaseDISCC/include/ReactionBaseDISCC.h
@@ -17,6 +17,8 @@
 // Define standard parameters used by SF and x-sections:
 #define BASE_PARS (int dataSetID, valarray<double> &val, map<string, valarray<double> > &err)
 
+class IntegrateDIS;
+
 class ReactionBaseDISCC : public ReactionTheory
 {
   public:
@@ -32,7 +34,7 @@ class ReactionBaseDISCC : public ReactionTheory
     virtual void setDatasetParameters( int dataSetID, map<string,string> pars, map<string,double> parsDataset) override ;
     virtual void initAtIteration() override;
     
-    virtual int compute(int dataSetID, valarray<double> &val, map<string, valarray<double> > &err);
+    virtual int compute(int dataSetID, valarray<double> &valExternal, map<string, valarray<double> > &errExternal);
   protected:
     enum class dataFlav { incl, c} ;      //!< Define final state.
 
@@ -72,6 +74,12 @@ class ReactionBaseDISCC : public ReactionTheory
     map <int,valarray<double> > _xf3u; 
     map <int,valarray<double> > _xf3d;
 
+  protected:
+    // for integrated cross sections
+    // method is based on legacy subroutine GetIntegratedDisXsection
+    map<int,IntegrateDIS*> _integrated;
+    virtual valarray<double> *GetBinValues(int idDS, const string& binName);
+
  protected:
     double _MW;
     double _Gf;
diff --git a/reactions/BaseDISCC/src/Makefile.am b/reactions/BaseDISCC/src/Makefile.am
index 8d717b9c0..bb2eea549 100644
--- a/reactions/BaseDISCC/src/Makefile.am
+++ b/reactions/BaseDISCC/src/Makefile.am
@@ -7,5 +7,7 @@ lib_LTLIBRARIES = libbasediscc_xfitter.la
 libbasediscc_xfitter_la_SOURCES = ReactionBaseDISCC.cc
 
 # libbasediscc_xfitter_la_LDFLAGS = place_if_needed  
+# link to main xFitter library (for IntegrateDIS)
+libbasediscc_xfitter_la_LDFLAGS = -lxfmain -L$(srcdir)/../../../src/.libs
 
 dist_noinst_HEADERS = ../include ../yaml
diff --git a/reactions/BaseDISCC/src/ReactionBaseDISCC.cc b/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
index 8430d02b9..6f1531573 100644
--- a/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
+++ b/reactions/BaseDISCC/src/ReactionBaseDISCC.cc
@@ -8,6 +8,7 @@
 
 #include "ReactionBaseDISCC.h"
 #include <iostream>
+#include <IntegrateDIS.h>
 
 // Helpers for QCDNUM (CC): 
 
@@ -15,8 +16,8 @@
 const double  CCEP2F[] = {0.,0.,1.,0.,1.,0., 0. ,1.,0.,1.,0.,0.,0.} ;
 const double  CCEM2F[] = {0.,0.,0.,1.,0.,1., 0. ,0.,1.,0.,1.,0.,0.} ;
 
-const double  CCEP3F[] = {0.,0.,-1.,0.,-1.,0.,0.,1.,0.,1.,0.,0.,0.};
-const double  CCEM3F[] = {0.,0. ,0.,-1.,0.,-1.,0.,0.,1.,0.,1.,0.,0.};
+const double  CCEP3F[] = {0.,0.,-1.,0.,-1.,0., 0., 1.,0.,1.,0.,0.,0.};
+const double  CCEM3F[] = {0.,0. ,0.,-1.,0.,-1., 0., 0.,1.,0.,1.,0.,0.};
 
 //! c
 // work in progress: according to 1001.2312 section 5,
@@ -65,8 +66,11 @@ int ReactionBaseDISCC::initAtStart(const string &s)
 }
 
 // Main function to compute results at an iteration
-int ReactionBaseDISCC::compute(int dataSetID, valarray<double> &val, map<string, valarray<double> > &err)
+int ReactionBaseDISCC::compute(int dataSetID, valarray<double> &valExternal, map<string, valarray<double> > &errExternal)
 {
+  valarray<double> val;
+  map<string, valarray<double> > err;
+
   // Basic formulat for CC cross section:
 
   auto *yp  = GetBinValues(dataSetID,"y");
@@ -84,6 +88,7 @@ int ReactionBaseDISCC::compute(int dataSetID, valarray<double> &val, map<string,
   FL (dataSetID,fl,err);
   xF3(dataSetID,xf3,err);
 
+
   double polarity = GetPolarisation(dataSetID);
   
   if ( GetCharge(dataSetID) > 0) {
@@ -95,6 +100,9 @@ int ReactionBaseDISCC::compute(int dataSetID, valarray<double> &val, map<string,
     val *= (1-polarity);
   }
 
+  //for(size_t i = 0; i < f2.size(); i++)
+  //  printf("%f %f    %f    %f    %f  =  %f\n", (*GetBinValues(dataSetID,"Q2"))[i], (*GetBinValues(dataSetID,"x"))[i], f2[i], fl[i], xf3[i], val[i]);
+
   if (! IsReduced(dataSetID)) {
     // extra factor for non-reduced cross section
     auto *xp  = GetBinValues(dataSetID,"x");
@@ -105,6 +113,20 @@ int ReactionBaseDISCC::compute(int dataSetID, valarray<double> &val, map<string,
     valarray<double> factor = (_MW*_MW*_MW*_MW/pow((q2+_MW*_MW),2))*_Gf*_Gf/(2*pi*x)*_convfac;
     val *= factor;
   }
+
+  if(_integrated.find(dataSetID) == _integrated.end())
+  {
+    // usual cross section at (q2,x) points
+    valExternal = val;
+    errExternal = err;
+  }
+  else
+  {
+    // integrated cross sections
+    valExternal = _integrated[dataSetID]->compute(val);
+    // no idea how error could be treated: for now do nothing
+    errExternal = err;
+  }
   
   return 0;
 }
@@ -127,33 +149,11 @@ void ReactionBaseDISCC::initAtIteration() {
 // 
 void  ReactionBaseDISCC::setDatasetParameters( int dataSetID, map<string,string> pars, map<string,double> parsDataset) 
 {
-  auto *q2p  = GetBinValues(dataSetID,"Q2"), *xp  = GetBinValues(dataSetID,"x"), *yp  = GetBinValues(dataSetID,"y");
-
-  // if Q2 and x bins and centre-of-mass energy provided, calculate y = Q2 / (s * x)
-  if(yp == nullptr && q2p != nullptr && xp != nullptr)
-  {
-    map<string,string>::iterator it = pars.find("energy");
-    if ( it != pars.end() )
-    {
-      double s = pow(stof(it->second), 2.0);
-      valarray<double> y = (*q2p) / (s * (*xp));
-      std::pair<string,valarray<double>* > dsBin = std::make_pair("y", &y);
-      AddBinning(dataSetID, &dsBin);
-      yp = GetBinValues(dataSetID, "y");
-    }
-  }
-
-  if (q2p == nullptr || xp == nullptr || yp == nullptr ) {
-    string msg = "F: Q2, x or Y bins are missing for CC DIS reaction for dataset " + std::to_string(dataSetID);
-    hf_errlog_(17100801,msg.c_str(), msg.size());
-  }
-  _npoints[dataSetID] = (*q2p).size();
   _polarisation[dataSetID] =  (parsDataset.find("epolarity") != parsDataset.end()) ? parsDataset["epolarity"] : 0;
   _charge[dataSetID]       =  (parsDataset.find("echarge")       != parsDataset.end()) ? parsDataset["echarge"] : 0;
   _isReduced[dataSetID]    =  (parsDataset.find("reduced")       != parsDataset.end()) ? parsDataset["reduced"] : 0;
 
   // check if settings are provided in the new format key=value
-
   // type: sigred, signonred (no F2, FL implemented so far, thus type is defined by bool _isReduced)
   // HERA data files provide 'signonred' CC cross sections
   // Inclusive "non-reduced" cross section by default.
@@ -161,7 +161,7 @@ void  ReactionBaseDISCC::setDatasetParameters( int dataSetID, map<string,string>
   string msg = "I: Calculating DIS CC reduced cross section";
   map<string,string>::iterator it = pars.find("type");
   if ( it != pars.end() ) {
-    if(it->second == "sigred") 
+    if(it->second == "sigred")
     {
       _isReduced[dataSetID] = 1;
       msg = "I: Calculating DIS CC reduced cross section";
@@ -209,7 +209,6 @@ void  ReactionBaseDISCC::setDatasetParameters( int dataSetID, map<string,string>
       hf_errlog_(18042502, str.c_str(), str.length());
     }
   }
-  hf_errlog_(17041001, msg.c_str(), msg.size());
 
   // e charge: double
   it = pars.find("echarge");
@@ -221,6 +220,65 @@ void  ReactionBaseDISCC::setDatasetParameters( int dataSetID, map<string,string>
   if ( it != pars.end() )
     _polarisation[dataSetID] = atof(it->second.c_str());
 
+  // check if centre-of-mass energy is provided
+  double s = -1.0;
+  map<string,string>::iterator itEnergy = pars.find("energy");
+  if ( itEnergy != pars.end() )
+    s = pow(stof(itEnergy->second), 2.0);
+
+  // bins
+  // if Q2min, Q2max, ymin and ymax (and optionally xmin, xmax) are provided, integrated cross sections are calculated
+  auto *q2minp  = GetBinValues(dataSetID,"Q2min");
+  auto *q2maxp  = GetBinValues(dataSetID,"Q2max");
+  // also try small first letter for Q2 (for backward compatibility)
+  if(!q2minp)
+    q2minp  = GetBinValues(dataSetID,"q2min");
+  if(!q2maxp)
+    q2maxp  = GetBinValues(dataSetID,"q2max");
+  auto *yminp  = GetBinValues(dataSetID,"ymin");
+  auto *ymaxp  = GetBinValues(dataSetID,"ymax");
+  // optional xmin, xmax for integrated cross sections
+  auto *xminp  = GetBinValues(dataSetID,"xmin");
+  auto *xmaxp  = GetBinValues(dataSetID,"xmax");
+
+  if(q2minp && q2maxp && yminp && ymaxp)
+  {
+    // integrated cross section
+    if(s < 0)
+      hf_errlog(18060100, "F: centre-of-mass energy is required for integrated DIS dataset " + std::to_string(dataSetID));
+    if(IsReduced(dataSetID))
+      hf_errlog(18060200, "F: integrated DIS can be calculated only for non-reduced cross sections, dataset " + std::to_string(dataSetID));
+    IntegrateDIS* iDIS = new IntegrateDIS();
+    _npoints[dataSetID] = iDIS->init(s, q2minp, q2maxp, yminp, ymaxp, xminp, xmaxp);
+    _integrated[dataSetID] = iDIS;
+    msg += " (integrated)";
+  }
+  else
+  {
+    // cross section at (Q2,x) points
+    auto *q2p  = GetBinValues(dataSetID,"Q2"), *xp  = GetBinValues(dataSetID,"x"), *yp  = GetBinValues(dataSetID,"y");
+
+    // if Q2 and x bins and centre-of-mass energy provided, calculate y = Q2 / (s * x)
+    if(yp == nullptr && q2p != nullptr && xp != nullptr)
+    {
+      if ( s > 0.0 )
+      {
+        valarray<double> y = (*q2p) / (s * (*xp));
+        std::pair<string,valarray<double>* > dsBin = std::make_pair("y", &y);
+        AddBinning(dataSetID, &dsBin);
+        yp = GetBinValues(dataSetID, "y");
+      }
+    }
+
+    if (q2p == nullptr || xp == nullptr || yp == nullptr ) {
+      string msg = "F: Q2, x or Y bins are missing for CC DIS reaction for dataset " + std::to_string(dataSetID);
+      hf_errlog_(17100801,msg.c_str(), msg.size());
+    }
+    _npoints[dataSetID] = (*q2p).size();
+  }
+
+  hf_errlog_(17041001, msg.c_str(), msg.size());
+
   // Allocate internal arrays:
   _f2u[dataSetID].resize(_npoints[dataSetID]);
   _f2d[dataSetID].resize(_npoints[dataSetID]);
@@ -230,6 +288,23 @@ void  ReactionBaseDISCC::setDatasetParameters( int dataSetID, map<string,string>
   _xf3d[dataSetID].resize(_npoints[dataSetID]);
 }
 
+valarray<double> *ReactionBaseDISCC::GetBinValues(int idDS, const string& binName)
+{
+  if(_integrated.find(idDS) == _integrated.end())
+    return ReactionTheory::GetBinValues(idDS, binName);
+  else
+  {
+    if(binName == "Q2")
+      return _integrated[idDS]->getBinValuesQ2();
+    else if(binName == "x")
+      return _integrated[idDS]->getBinValuesX();
+    else if(binName == "y")
+      return _integrated[idDS]->getBinValuesY();
+    else
+      return ReactionTheory::GetBinValues(idDS, binName);
+  }
+};
+
 
 // Get SF
 
@@ -291,6 +366,8 @@ void ReactionBaseDISCC::GetF2u(int dataSetID, valarray<double>& f2u)
         zmstfun_(id,CCEP2Fc[0], x[0], q2[0], (_f2u[dataSetID])[0], Npnt, flag);
         break ;
       }
+    //for(int i = 0; i < Npnt; i++)
+    //  printf("%f %f    %f\n", q2[i], x[i], (_f2u[dataSetID])[i]);
   }
   f2u = _f2u[dataSetID];
 }
@@ -335,10 +412,14 @@ void ReactionBaseDISCC::GetxF3u( int dataSetID, valarray<double>& xf3u )
         zmstfun_(id,CCEP3F[0], x[0], q2[0], (_xf3u[dataSetID])[0], Npnt, flag);
         break;
       case dataFlav::c :
+        //printf("before XF3u: %f\n", (_xf3u[dataSetID])[_xf3u[dataSetID].size() - 1]);
         zmstfun_(id,CCEP3Fc[0], x[0], q2[0], (_xf3u[dataSetID])[0], Npnt, flag);
+        //printf("after XF3u: %f\n", (_xf3u[dataSetID])[_xf3u[dataSetID].size() - 1]);
         break;
       }
+    //_xf3u[dataSetID] = _xf3u[dataSetID] * x;
   }
+  //printf("XF3u: %f\n", (_xf3u[dataSetID])[_xf3u[dataSetID].size() - 1]);
   xf3u = _xf3u[dataSetID];
 }
 
@@ -409,6 +490,7 @@ void ReactionBaseDISCC::GetxF3d( int dataSetID, valarray<double>& xf3d )
         zmstfun_(id,CCEM3Fc[0], x[0], q2[0], (_xf3d[dataSetID])[0], Npnt, flag);
         break;
     }
+    //_xf3d[dataSetID] = _xf3d[dataSetID] * x;
   }
   xf3d = _xf3d[dataSetID];
 }
diff --git a/reactions/BaseDISNC/include/ReactionBaseDISNC.h b/reactions/BaseDISNC/include/ReactionBaseDISNC.h
index b19136cc8..d6d455d73 100644
--- a/reactions/BaseDISNC/include/ReactionBaseDISNC.h
+++ b/reactions/BaseDISNC/include/ReactionBaseDISNC.h
@@ -16,6 +16,8 @@
 // Define standard parameters used by SF and x-sections:
 #define BASE_PARS (int dataSetID, valarray<double> &val, map<string, valarray<double> > &err)
 
+class IntegrateDIS;
+
 class ReactionBaseDISNC : public ReactionTheory
 {
  public:
@@ -29,7 +31,7 @@ class ReactionBaseDISNC : public ReactionTheory
     virtual void initAtIteration() override; 
     virtual int compute(int dataSetID, valarray<double> &val, map<string, valarray<double> > &err) override ;
  protected:
-    enum class dataType { sigred, f2, fl} ;  //!< Define compute output.
+    enum class dataType { signonred, sigred, f2, fl} ;  //!< Define compute output.
     enum class dataFlav { incl, c, b} ;      //!< Define final state.
 
     /* 
@@ -71,6 +73,7 @@ class ReactionBaseDISNC : public ReactionTheory
 
  protected:
     // some parameters which may change from iteration to iteration:
+    double _alphaem;
     double _Mz;
     double _Mw;
     double _sin2thetaW;
@@ -78,6 +81,9 @@ class ReactionBaseDISNC : public ReactionTheory
     double _au, _ad;
     double _vu, _vd;
 
+    // conversion constant factor
+    double _convfac;
+
  protected:
     const int GetNpoint(int dataSetID) {return _npoints[dataSetID];}
     const double GetPolarisation (int dataSetID) {return _polarisation[dataSetID];}
@@ -98,5 +104,11 @@ class ReactionBaseDISNC : public ReactionTheory
     map <int,valarray<double> > _fld; //!< FL for d-type quarks
     map <int,valarray<double> > _xf3u; 
     map <int,valarray<double> > _xf3d; 
+
+  protected:
+    // for integrated cross sections
+    // method is based on legacy subroutine GetIntegratedDisXsection
+    map<int,IntegrateDIS*> _integrated;
+    virtual valarray<double> *GetBinValues(int idDS, const string& binName);
 };
 
diff --git a/reactions/BaseDISNC/src/ReactionBaseDISNC.cc b/reactions/BaseDISNC/src/ReactionBaseDISNC.cc
index e26927300..a314d74e2 100644
--- a/reactions/BaseDISNC/src/ReactionBaseDISNC.cc
+++ b/reactions/BaseDISNC/src/ReactionBaseDISNC.cc
@@ -9,6 +9,7 @@
 #include "ReactionBaseDISNC.h"
 #include <iostream>
 #include <cstdio>
+#include <IntegrateDIS.h>
 
 template <typename T>
 void print(T d) {
@@ -47,14 +48,34 @@ extern "C" ReactionBaseDISNC* create() {
 // Initialize at the start of the computation
 int ReactionBaseDISNC::initAtStart(const string &s)
 {
+  _convfac = GetParam("convFac");
   return 0;
 }
 
 // Main function to compute results at an iteration
-int ReactionBaseDISNC::compute(int dataSetID, valarray<double> &val, map<string, valarray<double> > &err)
+int ReactionBaseDISNC::compute(int dataSetID, valarray<double> &valExternal, map<string, valarray<double> > &errExternal)
 {
+  valarray<double> val;
+  map<string, valarray<double> > err;
+
   switch ( GetDataType(dataSetID) )
     {
+    case dataType::signonred :
+      {
+        sred(dataSetID, val, err) ;
+        // transform reduced -> non-reduced cross sections
+        auto *xp  = GetBinValues(dataSetID,"x");
+        auto x = *xp;
+        auto *Q2p  = GetBinValues(dataSetID,"Q2");
+        auto q2 = *Q2p;
+        auto *yp  = GetBinValues(dataSetID,"y");
+        auto y = *yp;
+        const double pi = 3.1415926535897932384626433832795029;
+        valarray<double> yplus  = 1.0+(1.0-y)*(1.0-y);
+        valarray<double> factor = 2 * pi * _alphaem * _alphaem * yplus / (q2 * q2 * x) * _convfac;
+        val *= factor;
+        break ;
+      }
     case dataType::sigred :
       sred(dataSetID, val, err) ;
       break ;
@@ -65,10 +86,26 @@ int ReactionBaseDISNC::compute(int dataSetID, valarray<double> &val, map<string,
       FL(dataSetID, val, err) ;
       break ;
     }
+
+  if(_integrated.find(dataSetID) == _integrated.end())
+  {
+    // usual cross section at (q2,x) points
+    valExternal = val;
+    errExternal = err;
+  }
+  else
+  {
+    // integrated cross sections
+    valExternal = _integrated[dataSetID]->compute(val);
+    // no idea how error could be treated: for now do nothing
+    errExternal = err;
+  }
+
   return 0;
 }
 
 void ReactionBaseDISNC::initAtIteration() {
+  _alphaem = GetParam("alphaem");
   _Mz = GetParam("Mz");
   _Mw = GetParam("Mw");
   _sin2thetaW = GetParam("sin2thW");
@@ -93,27 +130,6 @@ void ReactionBaseDISNC::initAtIteration() {
 // 
 void  ReactionBaseDISNC::setDatasetParameters( int dataSetID, map<string,string> pars, map<string,double> parsDataset) 
 {
-  auto *q2p  = GetBinValues(dataSetID,"Q2"), *xp  = GetBinValues(dataSetID,"x"), *yp  = GetBinValues(dataSetID,"y");  
-
-  // if Q2 and x bins and centre-of-mass energy provided, calculate y = Q2 / (s * x)
-  if(yp == nullptr && q2p != nullptr && xp != nullptr)
-  {
-    map<string,string>::iterator it = pars.find("energy");
-    if ( it != pars.end() )
-    {
-      double s = pow(stof(it->second), 2.0);
-      valarray<double> y = (*q2p) / (s * (*xp));
-      std::pair<string,valarray<double>* > dsBin = std::make_pair("y", &y);
-      AddBinning(dataSetID, &dsBin);
-      yp = GetBinValues(dataSetID, "y");
-    }
-  }
-
-  if (q2p == nullptr || xp == nullptr || yp == nullptr ) {
-    string msg = "F: Q2, x or Y bins are missing for NC DIS reaction for dataset " + std::to_string(dataSetID);
-    hf_errlog_(17040801,msg.c_str(), msg.size());
-  }
-  _npoints[dataSetID] = (*q2p).size();
   _polarisation[dataSetID] =  (parsDataset.find("epolarity") != parsDataset.end()) ? parsDataset["epolarity"] : 0;
   _charge[dataSetID]       =  (parsDataset.find("echarge")       != parsDataset.end()) ? parsDataset["echarge"] : 0;
 
@@ -135,11 +151,15 @@ void  ReactionBaseDISNC::setDatasetParameters( int dataSetID, map<string,string>
   }
 
   // check if settings are provided in the new format key=value
-
-  // type: sigred, F2, FL
+  // type: signonred, sigred, F2, FL
   map<string,string>::iterator it = pars.find("type");
   if ( it != pars.end() ) {
-    if(it->second == "sigred")
+    if(it->second == "signonred")
+    {
+      _dataType[dataSetID] = dataType::signonred;
+      msg = "I: Calculating DIS NC double-differential (non-reduced) cross section";
+    }
+    else if(it->second == "sigred")
     {
       _dataType[dataSetID] = dataType::sigred;
       msg = "I: Calculating DIS NC reduced cross section";
@@ -189,7 +209,6 @@ void  ReactionBaseDISNC::setDatasetParameters( int dataSetID, map<string,string>
       hf_errlog_(17101902, str.c_str(), str.length());
     }
   }
-  hf_errlog_(17041001, msg.c_str(), msg.size());
 
   // e charge: double
   it = pars.find("echarge");
@@ -201,6 +220,65 @@ void  ReactionBaseDISNC::setDatasetParameters( int dataSetID, map<string,string>
   if ( it != pars.end() )
     _polarisation[dataSetID] = atof(it->second.c_str());
 
+  // check if centre-of-mass energy is provided
+  double s = -1.0;
+  map<string,string>::iterator itEnergy = pars.find("energy");
+  if ( itEnergy != pars.end() )
+    s = pow(stof(itEnergy->second), 2.0);
+
+  // bins
+  // if Q2min, Q2max, ymin and ymax (and optionally xmin, xmax) are provided, calculate integrated cross sections
+  auto *q2minp  = GetBinValues(dataSetID,"Q2min");
+  auto *q2maxp  = GetBinValues(dataSetID,"Q2max");
+  // also try small first letter for backward compatibility
+  if(!q2minp)
+    q2minp  = GetBinValues(dataSetID,"q2min");
+  if(!q2maxp)
+    q2maxp  = GetBinValues(dataSetID,"q2max");
+  auto *yminp  = GetBinValues(dataSetID,"ymin");
+  auto *ymaxp  = GetBinValues(dataSetID,"ymax");
+  // optional xmin, xmax for integrated cross sections
+  auto *xminp  = GetBinValues(dataSetID,"xmin");
+  auto *xmaxp  = GetBinValues(dataSetID,"xmax");
+
+  if(q2minp && q2maxp && yminp && ymaxp)
+  {
+    // integrated cross section
+    if(s < 0)
+      hf_errlog(18060100, "F: centre-of-mass energy is required for integrated DIS dataset " + std::to_string(dataSetID));
+    if(_dataType[dataSetID] != dataType::signonred)
+      hf_errlog(18060200, "F: integrated DIS can be calculated only for non-reduced cross sections, dataset " + std::to_string(dataSetID));
+    IntegrateDIS* iDIS = new IntegrateDIS();
+    _npoints[dataSetID] = iDIS->init(s, q2minp, q2maxp, yminp, ymaxp, xminp, xmaxp);
+    _integrated[dataSetID] = iDIS;
+    msg += " (integrated)";
+  }
+  else
+  {
+    // cross section at (Q2,x) points
+    auto *q2p  = GetBinValues(dataSetID,"Q2"), *xp  = GetBinValues(dataSetID,"x"), *yp  = GetBinValues(dataSetID,"y");
+
+    // if Q2 and x bins and centre-of-mass energy provided, calculate y = Q2 / (s * x)
+    if(yp == nullptr && q2p != nullptr && xp != nullptr)
+    {
+      if ( s > 0.0 )
+      {
+        valarray<double> y = (*q2p) / (s * (*xp));
+        std::pair<string,valarray<double>* > dsBin = std::make_pair("y", &y);
+        AddBinning(dataSetID, &dsBin);
+        yp = GetBinValues(dataSetID, "y");
+      }
+    }
+
+    if (q2p == nullptr || xp == nullptr || yp == nullptr ) {
+      string msg = "F: Q2, x or Y bins are missing for NC DIS reaction for dataset " + std::to_string(dataSetID);
+      hf_errlog_(17040801,msg.c_str(), msg.size());
+    }
+    _npoints[dataSetID] = (*q2p).size();
+  }
+
+  hf_errlog_(17041001, msg.c_str(), msg.size());
+
   // Allocate internal arrays:
   _f2u[dataSetID].resize(_npoints[dataSetID]);
   _f2d[dataSetID].resize(_npoints[dataSetID]);
@@ -210,6 +288,23 @@ void  ReactionBaseDISNC::setDatasetParameters( int dataSetID, map<string,string>
   _xf3d[dataSetID].resize(_npoints[dataSetID]);
 }
 
+valarray<double> *ReactionBaseDISNC::GetBinValues(int idDS, const string& binName)
+{
+  if(_integrated.find(idDS) == _integrated.end())
+    return ReactionTheory::GetBinValues(idDS, binName);
+  else
+  {
+    if(binName == "Q2")
+      return _integrated[idDS]->getBinValuesQ2();
+    else if(binName == "x")
+      return _integrated[idDS]->getBinValuesX();
+    else if(binName == "y")
+      return _integrated[idDS]->getBinValuesY();
+    else
+      return ReactionTheory::GetBinValues(idDS, binName);
+  }
+};
+
 void ReactionBaseDISNC::F2gamma BASE_PARS
 {
   valarray<double> f2u, f2d;
-- 
GitLab


From ab9e44704644d137f0589d536d360eee50d35f1e Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Sun, 3 Jun 2018 03:09:56 +0200
Subject: [PATCH 12/17] fix for non-reduced DIS NC: read convfac

---
 reactions/BaseDISNC/src/ReactionBaseDISNC.cc | 1 +
 1 file changed, 1 insertion(+)

diff --git a/reactions/BaseDISNC/src/ReactionBaseDISNC.cc b/reactions/BaseDISNC/src/ReactionBaseDISNC.cc
index a314d74e2..cd5b7e863 100644
--- a/reactions/BaseDISNC/src/ReactionBaseDISNC.cc
+++ b/reactions/BaseDISNC/src/ReactionBaseDISNC.cc
@@ -105,6 +105,7 @@ int ReactionBaseDISNC::compute(int dataSetID, valarray<double> &valExternal, map
 }
 
 void ReactionBaseDISNC::initAtIteration() {
+  _convfac = GetParam("convFac");
   _alphaem = GetParam("alphaem");
   _Mz = GetParam("Mz");
   _Mw = GetParam("Mw");
-- 
GitLab


From 163f2e2d1044506d0e95fee3f2630c61bcbc9c12 Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Sun, 3 Jun 2018 03:32:36 +0200
Subject: [PATCH 13/17] prevent division by y=0

---
 src/IntegrateDIS.cc | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/IntegrateDIS.cc b/src/IntegrateDIS.cc
index 8a2a1138e..9be77f87f 100644
--- a/src/IntegrateDIS.cc
+++ b/src/IntegrateDIS.cc
@@ -25,7 +25,9 @@ int IntegrateDIS::init(const double s,
     const double ymin = (*yminp)[i];
     const double ymax = (*ymaxp)[i];
 
-    //if(ymin == 0.)
+    // prevent 0 because there is division by y
+    if(ymin == 0.)
+      ymin = 1e-6;
 
     // just copy previous entry, if binning is the same
     //bool flagCopyValue = false;
-- 
GitLab


From 34d099703cc2d9f501790b3c08b0c051a71ad469 Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Sun, 3 Jun 2018 03:34:37 +0200
Subject: [PATCH 14/17] fixed compilation error

---
 src/IntegrateDIS.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/IntegrateDIS.cc b/src/IntegrateDIS.cc
index 9be77f87f..c29d8a4b9 100644
--- a/src/IntegrateDIS.cc
+++ b/src/IntegrateDIS.cc
@@ -22,7 +22,7 @@ int IntegrateDIS::init(const double s,
   {
     const double q2min = (*q2minp)[i];
     const double q2max = (*q2maxp)[i];
-    const double ymin = (*yminp)[i];
+    double ymin = (*yminp)[i];
     const double ymax = (*ymaxp)[i];
 
     // prevent 0 because there is division by y
-- 
GitLab


From 5fadaa974943203e04660e1dd21dc17c5e88ed3b Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Sun, 3 Jun 2018 18:11:16 +0200
Subject: [PATCH 15/17] call parent initAtIteration in FONLL (update parameters
 needed for non-reduced crosse sections

---
 reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc | 1 +
 reactions/FONLL_DISNC/src/ReactionFONLL_DISNC.cc | 1 +
 2 files changed, 2 insertions(+)

diff --git a/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc b/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
index b2a656ab6..47c0c37b3 100644
--- a/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
+++ b/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
@@ -136,6 +136,7 @@ int ReactionFONLL_DISCC::initAtStart(const string &s)
 // by the specific functions.
 void ReactionFONLL_DISCC::initAtIteration()
 {
+  ReactionBaseDISCC::initAtIteration ();
   // VB: With the following command, APFEL will be calling the "ExternalSetAPFEL1"
   // routine in FONLL/src/FONLL_wrap.f. This is not optimal but until that routine is
   // there, I cannot find a way to override it.
diff --git a/reactions/FONLL_DISNC/src/ReactionFONLL_DISNC.cc b/reactions/FONLL_DISNC/src/ReactionFONLL_DISNC.cc
index 0cdd38fd3..746b00d5f 100644
--- a/reactions/FONLL_DISNC/src/ReactionFONLL_DISNC.cc
+++ b/reactions/FONLL_DISNC/src/ReactionFONLL_DISNC.cc
@@ -131,6 +131,7 @@ int ReactionFONLL_DISNC::initAtStart(const string &s)
 // by the specific functions.
 void ReactionFONLL_DISNC::initAtIteration()
 {
+  ReactionBaseDISNC::initAtIteration ();
   // VB: With the following command, APFEL will be calling the "ExternalSetAPFEL1"
   // routine in FONLL/src/FONLL_wrap.f. This is not optimal but until that routine is
   // there, I cannot find a way to override it.
-- 
GitLab


From ae63ab2b50a0a64cfe92e903bf950c918c97ddda Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Mon, 2 Jul 2018 13:12:04 +0200
Subject: [PATCH 16/17] enabled FONLL at LO and with running HQ masses

---
 reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc | 14 +++++++++-----
 reactions/FONLL_DISNC/src/ReactionFONLL_DISNC.cc | 14 +++++++++-----
 2 files changed, 18 insertions(+), 10 deletions(-)

diff --git a/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc b/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
index 47c0c37b3..75f5ffe1e 100644
--- a/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
+++ b/reactions/FONLL_DISCC/src/ReactionFONLL_DISCC.cc
@@ -53,8 +53,8 @@ int ReactionFONLL_DISCC::initAtStart(const string &s)
 
   if (PtOrder == 0)
     {
-      const string msg = "F: FONLL at LO not available. Use the ZM-VFNS instead.";
-      hf_errlog_(17120601,msg.c_str(), msg.size());
+      //const string msg = "F: FONLL at LO not available. Use the ZM-VFNS instead.";
+      //hf_errlog_(17120601,msg.c_str(), msg.size());
     }
   else if (PtOrder == 1 && scheme == "FONLL-C")
     {
@@ -66,13 +66,18 @@ int ReactionFONLL_DISCC::initAtStart(const string &s)
       const string msg = "F: At NNLO only the FONLL-C scheme is possible";
       hf_errlog_(17120603,msg.c_str(), msg.size());
     }
+  else
+  {
+    APFEL::SetMassScheme(scheme);
+  }
 
   // If the MSbar masses are being used check that APFEL is used also
   // for the evolution.
   if (MassScheme == "MSbar")
     {
-      const string msg = "F: When using MSbar heavy quark masses it is necessarey to use APFEL for the DGLAP evolution";
-      hf_errlog_(17120604,msg.c_str(), msg.size());
+      // TODO check properly that APFEL is used for evolution
+      //const string msg = "F: When using MSbar heavy quark masses it is necessarey to use APFEL for the DGLAP evolution";
+      //hf_errlog_(17120604,msg.c_str(), msg.size());
     }
 
   // Set Parameters
@@ -93,7 +98,6 @@ int ReactionFONLL_DISCC::initAtStart(const string &s)
     }
   APFEL::SetAlphaQCDRef(Alphas_ref, Q_ref);
   APFEL::SetPerturbativeOrder(PtOrder);
-  APFEL::SetMassScheme(scheme);
   APFEL::SetNumberOfGrids(3);
   APFEL::SetGridParameters(1, 50, 3, 9.8e-7);
   APFEL::SetGridParameters(2, 40, 3, 1e-2);
diff --git a/reactions/FONLL_DISNC/src/ReactionFONLL_DISNC.cc b/reactions/FONLL_DISNC/src/ReactionFONLL_DISNC.cc
index 746b00d5f..02385f28f 100644
--- a/reactions/FONLL_DISNC/src/ReactionFONLL_DISNC.cc
+++ b/reactions/FONLL_DISNC/src/ReactionFONLL_DISNC.cc
@@ -48,8 +48,8 @@ int ReactionFONLL_DISNC::initAtStart(const string &s)
 
   if (PtOrder == 0)
     {
-      const string msg = "F: FONLL at LO not available. Use the ZM-VFNS instead.";
-      hf_errlog_(17120601,msg.c_str(), msg.size());
+      //const string msg = "F: FONLL at LO not available. Use the ZM-VFNS instead.";
+      //hf_errlog_(17120601,msg.c_str(), msg.size());
     }
   else if (PtOrder == 1 && scheme == "FONLL-C")
     {
@@ -61,13 +61,18 @@ int ReactionFONLL_DISNC::initAtStart(const string &s)
       const string msg = "F: At NNLO only the FONLL-C scheme is possible";
       hf_errlog_(17120603,msg.c_str(), msg.size());
     }
+  else
+  {
+    APFEL::SetMassScheme(scheme);
+  }
 
   // If the MSbar masses are being used check that APFEL is used also
   // for the evolution.
   if (MassScheme == "MSbar")
     {
-      const string msg = "F: When using MSbar heavy quark masses it is necessarey to use APFEL for the DGLAP evolution";
-      hf_errlog_(17120604,msg.c_str(), msg.size());
+      // TODO check properly that APFEL is used for evolution
+      //const string msg = "F: When using MSbar heavy quark masses it is necessarey to use APFEL for the DGLAP evolution";
+      //hf_errlog_(17120604,msg.c_str(), msg.size());
     }
 
   // Set Parameters
@@ -88,7 +93,6 @@ int ReactionFONLL_DISNC::initAtStart(const string &s)
     }
   APFEL::SetAlphaQCDRef(Alphas_ref, Q_ref);
   APFEL::SetPerturbativeOrder(PtOrder);
-  APFEL::SetMassScheme(scheme);
   APFEL::SetNumberOfGrids(3);
   APFEL::SetGridParameters(1, 50, 3, 9.8e-7);
   APFEL::SetGridParameters(2, 40, 3, 1e-2);
-- 
GitLab


From 2b780a95bb432b06b801d5caf6ad33360bb81cd2 Mon Sep 17 00:00:00 2001
From: Oleksandr Zenaiev <oleksandr.zenaiev@cern.ch>
Date: Wed, 18 Jul 2018 00:14:48 +0200
Subject: [PATCH 17/17] attempt to not link against lxfmain

---
 reactions/BaseDISCC/src/Makefile.am | 2 --
 1 file changed, 2 deletions(-)

diff --git a/reactions/BaseDISCC/src/Makefile.am b/reactions/BaseDISCC/src/Makefile.am
index bb2eea549..8d717b9c0 100644
--- a/reactions/BaseDISCC/src/Makefile.am
+++ b/reactions/BaseDISCC/src/Makefile.am
@@ -7,7 +7,5 @@ lib_LTLIBRARIES = libbasediscc_xfitter.la
 libbasediscc_xfitter_la_SOURCES = ReactionBaseDISCC.cc
 
 # libbasediscc_xfitter_la_LDFLAGS = place_if_needed  
-# link to main xFitter library (for IntegrateDIS)
-libbasediscc_xfitter_la_LDFLAGS = -lxfmain -L$(srcdir)/../../../src/.libs
 
 dist_noinst_HEADERS = ../include ../yaml
-- 
GitLab