diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
index 1995a7ea5d9208382bf1f21850c2317ce01bccfb..934c6513e131030854cc2c685850cc84c32ffb1e 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
@@ -372,7 +372,7 @@ bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol
double oxn1 = op->GetXNode1(obsBin, xNode + oOffsetXN1);
double rdiff = std::abs((txn1 - oxn1) / txn1);
if (rdiff > rtol) {
- debug["IsEquivalent"] << "rdiff too high: obsBin=" << obsBin << " xNode=" << xNode <<
+ debug["IsEquivalent"] << "XNode1 rdiff too high: obsBin=" << obsBin << " xNode=" << xNode <<
" txn1=" << txn1 << " oxn1=" << oxn1 << " rdiff=" << rdiff << endl;
return false;
}
@@ -393,33 +393,25 @@ bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol
for (unsigned int scaleNode = 0; scaleNode < tst3.size(); scaleNode++) {
fastNLO::v2d tst2 = tst3[scaleNode];
fastNLO::v2d ost2 = ost3[scaleNode];
- unsigned int tOffset2;
- unsigned int oOffset2;
- if (tst2.size() > ost2.size()) {
- tOffset2 = tst2.size() - ost2.size();
- oOffset2 = 0;
- } else {
- tOffset2 = 0;
- oOffset2 = ost2.size() - tst2.size();
- }
- for (unsigned int xNode = 0; xNode < std::min(tst2.size(), ost2.size()); xNode++) {
- fastNLO::v1d tst1 = tst2[xNode + tOffset2];
- fastNLO::v1d ost1 = ost2[xNode + oOffset2];
- if (tst1.size() != ost1.size()) {
- return false;
- }
- for (unsigned int subProcess = 0; subProcess < tst1.size(); subProcess++) {
- double rdiff = std::abs((tst1[subProcess] - ost1[subProcess]) / tst1[subProcess]);
- if (rdiff > rtol) {
- debug["IsEquivalent"] << "rdiff too high: obsBin=" << obsBin << " scaleVar=" << scaleVar
- << " scaleNode=" << scaleNode << " xNode=" << xNode << " subProcess=" << subProcess <<
- " t=" << tst1[subProcess] << " o=" << ost1[subProcess] << " rdiff=" << rdiff << endl;
+ for (int xHigh = 0; xHigh < std::min(tnb1, onb1); xHigh++) {
+ for (int xLow = 0; xLow <= xHigh; xLow++) {
+ fastNLO::v1d tst1 = tst2[GetXIndex(obsBin, xLow + tOffsetXN1, xHigh + tOffsetXN1)];
+ fastNLO::v1d ost1 = ost2[op->GetXIndex(obsBin, xLow + oOffsetXN1, xHigh + oOffsetXN1)];
+ if (tst1.size() != ost1.size()) {
return false;
}
+ for (unsigned int subProcess = 0; subProcess < tst1.size(); subProcess++) {
+ double rdiff = std::abs((tst1[subProcess] - ost1[subProcess]) / tst1[subProcess]);
+ if (rdiff > rtol) {
+ debug["IsEquivalent"] << "SigmaTilde rdiff too high: obsBin=" << obsBin << " scaleVar=" << scaleVar
+ << " scaleNode=" << scaleNode << " xHigh=" << xHigh << " xLow=" << xLow << " subProcess="
+ << subProcess << " t=" << tst1[subProcess] << " o=" << ost1[subProcess] << " rdiff="
+ << rdiff << endl;
+ return false;
+ }
+ }
}
-
}
-
}
}
}
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
index 1af1246c573cdd8bb35522edf7a0bfd959eeea33..f1dcbc9840687d76da9823b8867540ab43b70a5c 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
@@ -2470,7 +2470,12 @@ void fastNLOCreate::FillAllSubprocesses(const vector<vector<fnloEvent> >& events
double xmin = GetTheCoeffTable()->GetNPDFDim() == 1 ? std::min(fEvent._x1,fEvent._x2) : fEvent._x1;
double xmax = GetTheCoeffTable()->GetNPDFDim() == 1 ? std::max(fEvent._x1,fEvent._x2) : fEvent._x2;
vector<pair<int,double> > nxlo = fKernX1[ObsBin]->GetNodeValues(xmin);
- vector<pair<int,double> > nxup = fKernX2[ObsBin]->GetNodeValues(xmax);
+ vector<pair<int,double> > nxup;
+ if (c->NPDFDim > 1) {
+ nxup = fKernX2[ObsBin]->GetNodeValues(xmax);
+ } else {
+ nxup = fKernX1[ObsBin]->GetNodeValues(xmax);
+ }
fastNLO::v1d fg1 = fKernX1[ObsBin]->fgrid;
fastNLO::v1d xn1 = c->XNode1[ObsBin];
@@ -2492,19 +2497,66 @@ void fastNLOCreate::FillAllSubprocesses(const vector<vector<fnloEvent> >& events
fKernX1[ObsBin]->CheckX(xmin);
fKernX2[ObsBin]->CheckX(xmax);
ApplyPDFWeight(nxlo,xmin,fKernX1[ObsBin]->GetGridPtr()); // changes node values
- ApplyPDFWeight(nxup,xmax,fKernX2[ObsBin]->GetGridPtr()); // changes node values
+ if (c->NPDFDim > 1) {
+ ApplyPDFWeight(nxup,xmax,fKernX2[ObsBin]->GetGridPtr()); // changes node values
+ } else {
+ ApplyPDFWeight(nxup,xmax,fKernX1[ObsBin]->GetGridPtr()); // changes node values
+ }
}
+ // for (unsigned int i = 0; i < fg1.size(); i++) {
+ // std::cout << fg1[i] << " ";
+ // }
+ // std::cout << endl;
+ // for (unsigned int i = 0; i < nxlo.size(); i++) {
+ // std::cout << nxlo[i].first << " ";
+ // }
+ // std::cout << endl;
+ // for (unsigned int i = 0; i < nxup.size(); i++) {
+ // std::cout << fg1[nxlo[i].first] << " ";
+ // }
+ // std::cout << endl;
+ // for (unsigned int i = 0; i < nxup.size(); i++) {
+ // std::cout << nxlo[i].second << " ";
+ // }
+ // std::cout << endl;
+ // fastNLO::v1d fg2 = fKernX2[ObsBin]->fgrid;
+ // for (unsigned int i = 0; i < fg2.size(); i++) {
+ // std::cout << fg2[i] << " ";
+ // }
+ // std::cout << endl;
+ // for (unsigned int i = 0; i < nxup.size(); i++) {
+ // std::cout << nxup[i].first << " ";
+ // }
+ // std::cout << endl;
+ // for (unsigned int i = 0; i < nxup.size(); i++) {
+ // std::cout << fg2[nxup[i].first] << " ";
+ // }
+ // std::cout << endl;
+ // for (unsigned int i = 0; i < nxup.size(); i++) {
+ // std::cout << nxup[i].second << " ";
+ // }
+ // std::cout << endl;
fastNLO::v4d& st = c->SigmaTilde[ObsBin];
for (unsigned int is = 0 ; is<events.size() ; is++) {
double mu = fScenario._m1 * fScaleFac[is];
const vector<pair<int,double> >& nmu = fKernMuS[ObsBin][is]->GetNodeValues(mu);
for (unsigned int m1 = 0 ; m1<nmu.size() ; m1++) {
fastNLO::v2d& stm1 = st[is][nmu[m1].first];
- while(stm1.size() < maxNumNodes) {
- logger.debug["FillAllSubProcesses"] << "SigmaTilde ObsBin=" << ObsBin << " extended!" << endl;
- stm1.insert(stm1.begin(), fastNLO::v1d(NSubProc));
+ unsigned int i = 0;
+ while (stm1.size() < maxNumNodes) {
+ if (i * (i + 1) / 2 > stm1.size()) {
+ for (unsigned int j = 0; j < i; j++) {
+ stm1.insert(stm1.begin() + (j * (j + 1) / 2), fastNLO::v1d(NSubProc));
+ }
+ }
+ i++;
+ }
+ if (stm1.size() != maxNumNodes) {
+ logger.error["FillAllSubprocesses"] << "stm1 has size " << stm1.size()
+ << " when it should have size " << maxNumNodes << endl;
+ exit(1);
}
st[is][nmu[m1].first] = stm1;
for (unsigned int p = 0 ; p<events[is].size() ; p++) {
@@ -4532,7 +4584,10 @@ void fastNLOCreate::InitInterpolationKernels() {
// Create x grids with X_NNodes+1 nodes up to x_max = 1.
// The additional last node will be removed again below.
- int nxtot = fScenConsts.X_NNodes + 1;
+ int nxtot = fScenConsts.X_NNodes;
+ if (fScenConsts.X_NNodeCounting != "NodeDensity") {
+ nxtot += 1;
+ }
if (fScenConsts.X_NNodeCounting == "NodesPerMagnitude") { // "NodesMax","NodesPerBin","NodesPerMagnitude"
logger.debug["InitInterpolationKernels"]<<"Setting x nodes per magnitude: "<<fScenConsts.X_NNodeCounting<<endl;
fKernX1[i] = MakeInterpolationKernels(fScenConsts.X_Kernel,wrmX[i],1,fScenConsts.X_DistanceMeasure); // use 1 as upper x-value
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolBase.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolBase.cc
index c7a7bfce01b1e4ac7a8dfb218ada40af42db2128..8b1ccadc8baeadbe12a727c89237655675b0bb89 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolBase.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolBase.cc
@@ -30,9 +30,35 @@ fastNLOInterpolBase::fastNLOInterpolBase(double density, fastNLOGrid::GridType t
fLastGridPointWasRemoved=false;
debug["fastNLOInterpolBase"]<<"Distance measure = "<<type<<endl;
fdm = type;
+ fExtendLow = true;
fvalmax = 1.0;
- fvalmin = Function_sqrtlog10_inv(-nMinNodes/density);
+ double valMinH = -(nMinNodes - 1) / density;
+ switch (fdm) {
+ case fastNLOGrid::kLinear:
+ fvalmin = valMinH;
+ break;
+ case fastNLOGrid::kLogLog025:
+ fvalmin = Function_loglog025_inv(valMinH);
+ break;
+ case fastNLOGrid::kLog10:
+ fvalmin = Function_log10_inv(valMinH);
+ break;
+ case fastNLOGrid::kSqrtLog10:
+ fvalmin = Function_sqrtlog10_inv(valMinH);
+ break;
+ case fastNLOGrid::kLogLog:
+ fvalmin = Function_loglog_inv(valMinH);
+ break;
+ case fastNLOGrid::k3rdrtLog10:
+ fvalmin = Function_3rdrtlog10_inv(valMinH);
+ break;
+ case fastNLOGrid::k4thrtLog10:
+ fvalmin = Function_4thrtlog10_inv(valMinH);
+ break;
+ default:
+ error["MakeGridFromHGrid"]<<"Unknown grid type."<<endl;
+ }
}
@@ -141,24 +167,28 @@ void fastNLOInterpolBase::MakeGrids(int nNodes, double ReduceXmin){
}
-int fastNLOInterpolBase::FindLargestPossibleNode(double x){
+int fastNLOInterpolBase::FindLargestPossibleNode(double x, bool canExtend = false){
// --- find x position in range: 0 <= node1 < nnode1-1
int node1 = fgrid.size()-2; // --- initialize with largest possible value
if ( fLastGridPointWasRemoved ) node1=fgrid.size()-1;
- bool gridExtended = false;
- while ( x < fgrid[0] ) {
- double newxH = fHgrid[0] - (fHgrid[1] - fHgrid[0]);
- fHgrid.insert(fHgrid.begin(), newxH);
- fgrid.insert(fgrid.begin(), Function_sqrtlog10_inv(newxH));
- debug["FindLargestPossibleNode"]<<"Value is smaller than smallest node. Extending grid. x="<<x<<endl;
- //warn["FindLargestPossibleNode"]<<"Value is smaller than smallest node. Using first node. This may bias the result! x="<<x<<endl;
- gridExtended = true;
- }
- if (gridExtended) {
- return 0;
+ if (canExtend && fExtendLow) {
+ bool gridExtended = false;
+ while ( x < fgrid[1] ) {
+ double newxH = fHgrid[0] - (fHgrid[1] - fHgrid[0]);
+ fHgrid.insert(fHgrid.begin(), newxH);
+ fgrid.insert(fgrid.begin(), Function_sqrtlog10_inv(newxH));
+ debug["FindLargestPossibleNode"]<<"Value is smaller than smallest node. Extending grid. x="<<x<<endl;
+ gridExtended = true;
+ }
+ if (gridExtended) {
+ return 1;
+ }
}
if ( x==fgrid[0] ) {
return 0;
+ } else if ( x < fgrid[0] ) {
+ warn["FindLargestPossibleNode"]<<"Value is smaller than smallest node. Using first node. This may bias the result! x="<<x<<endl;
+ return 0;
}
if ( x > fgrid.back() ) {
if ( !fLastGridPointWasRemoved )
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolCatmullRom.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolCatmullRom.cc
index 45536a92fd18e8c1751af0f1d8965840b8691e87..db59a7acebfc7181b1b441296883cf96202fad9b 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolCatmullRom.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolCatmullRom.cc
@@ -52,7 +52,7 @@ void fastNLOInterpolCatmullRom::CalcNodeValues(vector<pair<int,double> >& nodes,
// --- get scale interpolation kernel and updated scalenode position: 1 <= nmu < ntot-2
int nmod = 0; // --- variable for final node
- int nnode = FindLargestPossibleNode(x);
+ int nnode = FindLargestPossibleNode(x, true);
int nmax = fgrid.size()-2;
if (fLastGridPointWasRemoved) nmax=fgrid.size()-1;
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolLagrange.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolLagrange.cc
index e9661eec1923c0e4dfb1f1c66d28b3db67afe292..5458be9bbab3ce55463fa947746dfc292f6baff5 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolLagrange.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolLagrange.cc
@@ -39,7 +39,7 @@ void fastNLOInterpolLagrange::CalcNodeValues(vector<pair<int,double> >& nodes, d
static const unsigned int nS = 4; // number of nodes that receive contributions from this interpolation
- int nnode = FindLargestPossibleNode(x);
+ int nnode = FindLargestPossibleNode(x, true);
// --- relative distance delta - in function fdm H(x)
// deltascale (Interpol(.,.,.delta,.): relative distance of value to node 'nnode'
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolLinear.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolLinear.cc
index fceef1212752f4f9e5c793798d7b971c0856f057..0ea9c793e4c43e388e39251363535c9931a035ca 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolLinear.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolLinear.cc
@@ -35,7 +35,7 @@ void fastNLOInterpolLinear::CalcNodeValues(vector<pair<int,double> >& nodes, dou
// deltascale (Interpol(.,.,.delta,.): relative distance of value to node 'nnode'
double delta = GetDelta(x);
// --- get scale interpolation kernel and updated scalenode position: 1 <= nmu < ntot-2
- int nnode = FindLargestPossibleNode(x);
+ int nnode = FindLargestPossibleNode(x, true);
// --- set nodes
nodes.resize(2);
nodes[0] = make_pair(nnode, 1.-delta);
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOInterpolBase.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOInterpolBase.h
index b019c96058fdd2d397c2d09db7ef40d127ab101f..9a7169e3fdd63f0a74919c112ccdc690b808b3fe 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOInterpolBase.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOInterpolBase.h
@@ -59,7 +59,7 @@ protected:
//virtual std::vector<std::pair<int,double> > CalcNodeValues(double val) = 0;
virtual void CalcNodeValues(std::vector<std::pair<int,double> >& nodes, double val) = 0;
- int FindLargestPossibleNode(double);
+ int FindLargestPossibleNode(double, bool);
inline double Function_loglog025( double mu ){
// function H(mu) = log(log( mu / 0.25 ))
@@ -108,6 +108,7 @@ protected:
fastNLOGrid::GridType fdm; // distance measure
std::vector<double> fHgrid;
int fnmod ; // variable for final nodes. Has to be filled by inherited algorithm
+ bool fExtendLow = false;
};