From c78f89c883512c7ebdee93ec76716d3d6a950f36 Mon Sep 17 00:00:00 2001
From: JohannesGaessler <johannesg@5d6.de>
Date: Sun, 26 Feb 2023 21:55:44 +0100
Subject: [PATCH] Bin density implementation for flexible scale
---
.../fastnlotoolkit/fastNLOCoeffAddFix.cc | 10 +-
.../fastnlotoolkit/fastNLOCoeffAddFlex.cc | 102 ++++++++----------
v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc | 56 ++++++----
v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc | 10 +-
.../include/fastnlotk/fastNLOCoeffAddFlex.h | 2 +-
.../include/fastnlotk/fastNLOTools.h | 3 +-
6 files changed, 96 insertions(+), 87 deletions(-)
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
index 964a6a0a..db2d620f 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
@@ -362,6 +362,7 @@ bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol
}
fastNLO::v5d ost5 = op->SigmaTilde;
if (SigmaTilde.size() != ost5.size()) {
+ debug["IsEquivalent"] << "Number of observable bins is different." << endl;
return false;
}
for (unsigned int obsBin = 0; obsBin < SigmaTilde.size(); obsBin++) {
@@ -369,6 +370,7 @@ bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol
unsigned int oOffsetXN1;
int tnb1 = GetNxtot1(obsBin);
int onb1 = op->GetNxtot1(obsBin);
+ int xMax = std::min(tnb1, onb1);
if (tnb1 > onb1) {
tOffsetXN1 = tnb1 - onb1;
oOffsetXN1 = 0;
@@ -391,7 +393,7 @@ 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];
- for (int xHigh = 0; xHigh < std::min(tnb1, onb1); xHigh++) {
+ for (int xHigh = 0; xHigh < xMax; 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)];
@@ -402,9 +404,9 @@ bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol
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;
+ << " scaleNode=" << scaleNode << " xHigh=" << (xHigh - xMax) << " xLow=" << (xLow - xMax)
+ << " subProcess=" << subProcess << " t=" << tst1[subProcess] << " o=" << ost1[subProcess]
+ << " rdiff=" << rdiff << endl;
return false;
}
}
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFlex.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFlex.cc
index a8b03397..47637be9 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFlex.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFlex.cc
@@ -583,31 +583,17 @@ void fastNLOCoeffAddFlex::CatBin(const fastNLOCoeffAddFlex& other, unsigned int
bool fastNLOCoeffAddFlex::IsEquivalent(const fastNLOCoeffBase& other, double rtol) const {
const fastNLOCoeffAddFlex* op = dynamic_cast<const fastNLOCoeffAddFlex*>(&other);
if (op == nullptr) {
+ debug["IsEquivalent"] << "other is not of type fastNLOCoeffAddFlex." << endl;
return false;
}
- for (unsigned int obsBin = 0; obsBin < SigmaTildeMuIndep.size(); obsBin++) {
- unsigned int tOffsetXN1;
- unsigned int oOffsetXN1;
- int tnb1 = GetNxtot1(obsBin);
- int onb1 = op->GetNxtot1(obsBin);
- if (tnb1 > onb1) {
- tOffsetXN1 = tnb1 - onb1;
- oOffsetXN1 = 0;
- } else {
- tOffsetXN1 = 0;
- oOffsetXN1 = onb1 - tnb1;
- }
- for (int xNode = 0; xNode < std::min(tnb1, onb1); xNode++) {
- double txn1 = GetXNode1(obsBin, xNode + tOffsetXN1);
- 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 <<
- " txn1=" << txn1 << " oxn1=" << oxn1 << " rdiff=" << rdiff << endl;
- return false;
- }
- }
+ if (!fastNLOTools::SameTails(GetAllXNodes1(), op->GetAllXNodes1(), rtol)) {
+ debug["IsEquivalent"] << "XNode1 not equivalent, see above." << endl;
+ return false;
+ }
+ if (!fastNLOTools::SameTails(GetAllXNodes2(), op->GetAllXNodes2(), rtol)) {
+ debug["IsEquivalent"] << "XNode2 not equivalent, see above." << endl;
+ return false;
}
std::vector<const fastNLO::v5d*> tst6 = GetSigmaTildes();
std::vector<const fastNLO::v5d*> ost6 = op->GetSigmaTildes();
@@ -615,7 +601,7 @@ bool fastNLOCoeffAddFlex::IsEquivalent(const fastNLOCoeffBase& other, double rto
for (unsigned int i = 0; i < tst6.size(); i++) {
const fastNLO::v5d* tst5 = tst6[i];
const fastNLO::v5d* ost5 = ost6[i];
- if (!IsSigmaTildeEquivalent(tst5, ost5, rtol, names[i])) {
+ if (!IsSigmaTildeEquivalent(op, tst5, ost5, rtol, names[i])) {
return false;
}
}
@@ -624,53 +610,57 @@ bool fastNLOCoeffAddFlex::IsEquivalent(const fastNLOCoeffBase& other, double rto
//________________________________________________________________________________________________________________ //
-bool fastNLOCoeffAddFlex::IsSigmaTildeEquivalent(const fastNLO::v5d *tst5, const fastNLO::v5d *ost5, double rtol, std::string name) const {
+bool fastNLOCoeffAddFlex::IsSigmaTildeEquivalent(const fastNLOCoeffAddFlex* op, const fastNLO::v5d *tst5, const fastNLO::v5d *ost5, double rtol, std::string name) const {
if (tst5->size() != ost5->size()) {
+ debug["IsEquivalent"] << name << ": Number of observable bins is different." << endl;
return false;
}
for (unsigned int obsBin = 0; obsBin < tst5->size(); obsBin++) {
-
fastNLO::v4d tst4 = tst5->at(obsBin);
fastNLO::v4d ost4 = ost5->at(obsBin);
- if (tst4.size() != ost4.size()) {
- return false;
+ unsigned int tOffsetXN1;
+ unsigned int oOffsetXN1;
+ int tnb1 = GetNxtot1(obsBin);
+ int onb1 = op->GetNxtot1(obsBin);
+ int xMax = std::min(tnb1, onb1);
+ if (tnb1 > onb1) {
+ tOffsetXN1 = tnb1 - onb1;
+ oOffsetXN1 = 0;
+ } else {
+ tOffsetXN1 = 0;
+ oOffsetXN1 = onb1 - tnb1;
}
- for (unsigned int scaleVar = 0; scaleVar < tst4.size(); scaleVar++) {
- fastNLO::v3d tst3 = tst4[scaleVar];
- fastNLO::v3d ost3 = ost4[scaleVar];
- if (tst3.size() != ost3.size()) {
- return false;
- }
- 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 (int xHigh = 0; xHigh < xMax; xHigh++) {
+ for (int xLow = 0; xLow <= xHigh; xLow++) {
+ fastNLO::v3d tst3 = tst4[GetXIndex(obsBin, xLow + tOffsetXN1, xHigh + tOffsetXN1)];
+ fastNLO::v3d ost3 = ost4[op->GetXIndex(obsBin, xLow + oOffsetXN1, xHigh + oOffsetXN1)];
+ if (tst3.size() != ost3.size()) {
+ return false;
}
- 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()) {
+ for (unsigned int scaleNode1 = 0; scaleNode1 < tst3.size(); scaleNode1++) {
+ fastNLO::v2d tst2 = tst3[scaleNode1];
+ fastNLO::v2d ost2 = ost3[scaleNode1];
+ if (tst2.size() != ost2.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"] << name << ": rdiff too high: obsBin=" << obsBin << " scaleVar=" << scaleVar
- << " scaleNode=" << scaleNode << " xNode=" << xNode << " subProcess=" << subProcess <<
- " t=" << tst1[subProcess] << " o=" << ost1[subProcess] << " rdiff=" << rdiff << endl;
+ for (unsigned int scaleNode2 = 0; scaleNode2 < tst2.size(); scaleNode2++) {
+ fastNLO::v1d tst1 = tst2[scaleNode2];
+ fastNLO::v1d ost1 = ost2[scaleNode2];
+ 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 << " scaleNode1=" << scaleNode1
+ << " scaleNode2=" << scaleNode2 << " xHigh=" << (xHigh - xMax) << " xLow=" << (xLow - xMax)
+ << " 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 fc124a70..885e3c83 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
@@ -2497,7 +2497,8 @@ void fastNLOCreate::FillAllSubprocesses(const vector<vector<fnloEvent> >& events
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];
- fastNLOTools::ExtendSigmaTildeX(stm1, c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, NSubProc);
+ fastNLOTools::ExtendSigmaTildeX(
+ stm1, c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, fastNLO::v1d(NSubProc));
for (unsigned int p = 0 ; p<events[is].size() ; p++) {
double wgt = wgtfac * events[is][p]._w * nmu[m1].second / BinSize[ObsBin];
// .......................................................................................
@@ -3004,8 +3005,16 @@ void fastNLOCreate::FillContributionFlexHHC(fastNLOCoeffAddFlex* c, int ObsBin,
//cout<<"try to interpol. ObsBin="<<ObsBin<<" ,x1="<<fEvent._x1<<", x2="<<fEvent._x2<<", mu1="<<Scenario._m1<<", mu2="<<Scenario._m2<<endl;
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> > nxup;
vector<pair<int,double> > nxlo = fKernX1[ObsBin]->GetNodeValues(xmin);
- vector<pair<int,double> > nxup = fKernX2[ObsBin]->GetNodeValues(xmax);
+ if (c->NPDFDim > 1) {
+ nxup = fKernX2[ObsBin]->GetNodeValues(xmax);
+ fastNLOTools::ExtendHead(c->XNode1[ObsBin], fKernX1[ObsBin]->fgrid);
+ fastNLOTools::ExtendHead(c->XNode2[ObsBin], fKernX2[ObsBin]->fgrid);
+ } else {
+ nxup = fKernX1[ObsBin]->GetNodeValues(xmax);
+ fastNLOTools::ExtendHead(c->XNode1[ObsBin], fKernX1[ObsBin]->fgrid);
+ }
const vector<pair<int,double> >& nmu1 = fKernMu1[ObsBin]->GetNodeValues(fScenario._m1);
const vector<pair<int,double> >& nmu2 = fKernMu2[ObsBin]->GetNodeValues(fScenario._m2);
@@ -3013,7 +3022,11 @@ void fastNLOCreate::FillContributionFlexHHC(fastNLOCoeffAddFlex* c, int ObsBin,
fKernX1[ObsBin]->CheckX(xmin);
fKernX2[ObsBin]->CheckX(xmax);
ApplyPDFWeight(nxlo,xmin,fKernX1[ObsBin]->GetGridPtr());
- ApplyPDFWeight(nxup,xmax,fKernX2[ObsBin]->GetGridPtr());
+ if (c->NPDFDim > 1) {
+ ApplyPDFWeight(nxup,xmax,fKernX2[ObsBin]->GetGridPtr());
+ } else {
+ ApplyPDFWeight(nxup,xmax,fKernX1[ObsBin]->GetGridPtr());
+ }
}
@@ -3027,20 +3040,6 @@ void fastNLOCreate::FillContributionFlexHHC(fastNLOCoeffAddFlex* c, int ObsBin,
HalfMatrixCheck(fEvent._x1,fEvent._x2,xminbin,xmaxbin,p);
//HalfMatrixCheck(xminbin,xmaxbin,p);
int ixHM = GetXIndex(ObsBin,xminbin,xmaxbin);
- fastNLO::v1d fg1 = fKernX1[ObsBin]->fgrid;
- fastNLO::v1d xn1 = c->XNode1[ObsBin];
- for (unsigned int i = fg1.size() - xn1.size(); i > 0; i--) {
- xn1.insert(xn1.begin(), fg1[i - 1]);
- }
- c->XNode1[ObsBin] = xn1;
- if (c->NPDFDim > 1) {
- fastNLO::v1d fg2 = fKernX2[ObsBin]->fgrid;
- fastNLO::v1d xn2 = c->XNode2[ObsBin];
- for (unsigned int i = fg2.size() - xn2.size(); i > 0; i--) {
- xn2.insert(xn2.begin(), fg2[i - 1]);
- }
- c->XNode2[ObsBin] = xn2;
- }
ExtendAllFlexSigmaTilde(c, ObsBin);
for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) {
@@ -3243,12 +3242,23 @@ void fastNLOCreate::FillContributionFlexDIS(fastNLOCoeffAddFlex* c, int ObsBin,
// ___________________________________________________________________________________________________
void fastNLOCreate::ExtendAllFlexSigmaTilde(fastNLOCoeffAddFlex* c, int ObsBin) {
- ExtendFlexSigmaTilde(c, ObsBin, c->SigmaTildeMuIndep);
- ExtendFlexSigmaTilde(c, ObsBin, c->SigmaTildeMuFDep);
- ExtendFlexSigmaTilde(c, ObsBin, c->SigmaTildeMuRDep);
- ExtendFlexSigmaTilde(c, ObsBin, c->SigmaTildeMuRRDep);
- ExtendFlexSigmaTilde(c, ObsBin, c->SigmaTildeMuFFDep);
- ExtendFlexSigmaTilde(c, ObsBin, c->SigmaTildeMuRFDep);
+ fastNLO::v3d insertValue = fastNLO::v3d(
+ c->GetNScaleNode1(ObsBin),
+ fastNLO::v2d(
+ c->GetNScaleNode1(ObsBin),
+ fastNLO::v1d(c->GetNSubproc())));
+ fastNLOTools::ExtendSigmaTildeX(
+ c->SigmaTildeMuIndep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+ fastNLOTools::ExtendSigmaTildeX(
+ c->SigmaTildeMuRDep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+ fastNLOTools::ExtendSigmaTildeX(
+ c->SigmaTildeMuFDep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+ fastNLOTools::ExtendSigmaTildeX(
+ c->SigmaTildeMuRRDep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+ fastNLOTools::ExtendSigmaTildeX(
+ c->SigmaTildeMuFFDep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+ fastNLOTools::ExtendSigmaTildeX(
+ c->SigmaTildeMuRFDep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
}
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc
index 2093f9cf..973b7fdf 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc
@@ -580,13 +580,17 @@ namespace fastNLOTools {
}
//______________________________________________________________________________
- void ExtendSigmaTildeX(fastNLO::v2d& SigmaTildeX, unsigned int DimSize1, unsigned int DimSize2, int NPDFDim, int NSubProc) {
+ template <typename T> void ExtendSigmaTildeX(
+ std::vector<T>& SigmaTildeX, unsigned int DimSize1, unsigned int DimSize2, int NPDFDim, T InsertValue) {
+ if (SigmaTildeX.size() >= DimSize1 * (DimSize1 + 1) / 2) {
+ return;
+ }
if (NPDFDim == 1) {
unsigned int i = 1;
while (SigmaTildeX.size() < DimSize1 * (DimSize1 + 1) / 2) {
if (i * (i + 1) / 2 > SigmaTildeX.size()) {
for (unsigned int j = 0; j < i; j++) {
- SigmaTildeX.insert(SigmaTildeX.begin() + (j * (j + 1) / 2), fastNLO::v1d(NSubProc));
+ SigmaTildeX.insert(SigmaTildeX.begin() + (j * (j + 1) / 2), InsertValue);
}
}
i++;
@@ -596,5 +600,7 @@ namespace fastNLOTools {
exit(1);
}
}
+ template void ExtendSigmaTildeX<fastNLO::v1d>(fastNLO::v2d&, unsigned int, unsigned int, int, fastNLO::v1d);
+ template void ExtendSigmaTildeX<fastNLO::v3d>(fastNLO::v4d&, unsigned int, unsigned int, int, fastNLO::v3d);
} // end namespace fastNLO
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddFlex.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddFlex.h
index 9326c03c..d3b2e88d 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddFlex.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddFlex.h
@@ -49,7 +49,7 @@ public:
return {&SigmaTildeMuIndep,&SigmaTildeMuRDep,&SigmaTildeMuFDep,&SigmaTildeMuRRDep,&SigmaTildeMuFFDep,&SigmaTildeMuRFDep};
};//!< Get access to sigma tilde
bool IsEquivalent(const fastNLOCoeffBase& other, double rtol) const;
- bool IsSigmaTildeEquivalent(const fastNLO::v5d *tst5, const fastNLO::v5d *ost5, double rtol, std::string name) const;
+ bool IsSigmaTildeEquivalent(const fastNLOCoeffAddFlex* op, const fastNLO::v5d *tst5, const fastNLO::v5d *ost5, double rtol, std::string name) const;
protected:
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h
index a3d3f548..b5c7f693 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h
@@ -94,7 +94,8 @@ namespace fastNLOTools {
bool SameTails(fastNLO::v1d vector1, fastNLO::v1d vector2, double rtol=0.0);
bool SameTails(fastNLO::v2d vector1, fastNLO::v2d vector2, double rtol=0.0);
void ExtendHead(fastNLO::v1d& vector1, fastNLO::v1d& vector2);
- void ExtendSigmaTildeX(fastNLO::v2d& SigmaTildeX, unsigned int DimSize1, unsigned int DimSize2, int NPDFDim, int NSubProc);
+ template <typename T> void ExtendSigmaTildeX(
+ std::vector<T>& SigmaTildeX, unsigned int DimSize1, unsigned int DimSize2, int NPDFDim, T InsertValue);
};
--
GitLab