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