diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc index 6947f5fdd796aa4dd81c0178f29aa345166a79b0..82d490bb8db8f07bb02a98c175ab10b605af7a6d 100644 --- a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc +++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc @@ -2466,19 +2466,37 @@ void fastNLOCreate::FillAllSubprocesses(const vector<vector<fnloEvent> >& events fastNLOCoeffAddFix* c = (fastNLOCoeffAddFix*)GetTheCoeffTable(); int NSubProc = c->GetNSubproc(); + int NPDFDim = c->GetNPDFDim(); + int NScaleVar = c->GetNScalevar(); + int NScaleNode = c->GetNScaleNode(); // do interpolation - double xmin = c->GetNPDFDim() == 1 ? std::min(fEvent._x1,fEvent._x2) : fEvent._x1; - double xmax = c->GetNPDFDim() == 1 ? std::max(fEvent._x1,fEvent._x2) : fEvent._x2; + double xmin = NPDFDim == 1 ? std::min(fEvent._x1,fEvent._x2) : fEvent._x1; + double xmax = NPDFDim == 1 ? std::max(fEvent._x1,fEvent._x2) : fEvent._x2; vector<pair<int,double> > nxlo = fKernX1[ObsBin]->GetNodeValues(xmin); vector<pair<int,double> > nxup; - if (c->NPDFDim > 1) { + unsigned int oldNxtot1 = c->XNode1[ObsBin].size(); + unsigned int oldNxtot2, newNxtot2; + if (NPDFDim > 1) { nxup = fKernX2[ObsBin]->GetNodeValues(xmax); fastNLOTools::ExtendHead(c->XNode1[ObsBin], fKernX1[ObsBin]->fgrid); + oldNxtot2 = c->XNode2[ObsBin].size(); fastNLOTools::ExtendHead(c->XNode2[ObsBin], fKernX2[ObsBin]->fgrid); + newNxtot2 = c->XNode2[ObsBin].size(); } else { nxup = fKernX1[ObsBin]->GetNodeValues(xmax); fastNLOTools::ExtendHead(c->XNode1[ObsBin], fKernX1[ObsBin]->fgrid); } + unsigned int newNxtot1 = c->XNode1[ObsBin].size(); + fastNLO::v4d& st = c->SigmaTilde[ObsBin]; + if (newNxtot1 > oldNxtot1 || (NPDFDim > 1 && newNxtot2 > oldNxtot2)) { + for (int scaleVar = 0; scaleVar < NScaleVar; scaleVar++) { + for (int scaleNode = 0; scaleNode < NScaleNode; scaleNode++) { + fastNLOTools::ExtendSigmaTildeX( + st[scaleVar][scaleNode], oldNxtot1, newNxtot1, oldNxtot2, newNxtot2, + NPDFDim, fastNLO::v1d(NSubProc)); + } + } + } if (fApplyPDFReweight) { fKernX1[ObsBin]->CheckX(xmin); @@ -2491,14 +2509,11 @@ void fastNLOCreate::FillAllSubprocesses(const vector<vector<fnloEvent> >& events } } - 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]; - 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]; // ....................................................................................... @@ -3007,6 +3022,8 @@ void fastNLOCreate::FillContributionFlexHHC(fastNLOCoeffAddFlex* c, int ObsBin, 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); + unsigned int oldNxtot1 = c->GetNxtot1(ObsBin); + unsigned int oldNxtot2 = c->GetNxtot2(ObsBin); if (c->NPDFDim > 1) { nxup = fKernX2[ObsBin]->GetNodeValues(xmax); fastNLOTools::ExtendHead(c->XNode1[ObsBin], fKernX1[ObsBin]->fgrid); @@ -3040,7 +3057,7 @@ 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); - ExtendAllFlexSigmaTilde(c, ObsBin); + ExtendAllFlexSigmaTilde(c, oldNxtot1, oldNxtot2, ObsBin); for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) { for (unsigned int mu2 = 0 ; mu2<nmu2.size() ; mu2++) { @@ -3241,24 +3258,31 @@ void fastNLOCreate::FillContributionFlexDIS(fastNLOCoeffAddFlex* c, int ObsBin, // ___________________________________________________________________________________________________ -void fastNLOCreate::ExtendAllFlexSigmaTilde(fastNLOCoeffAddFlex* c, int ObsBin) { +void fastNLOCreate::ExtendAllFlexSigmaTilde( + fastNLOCoeffAddFlex* c, unsigned int OldNxtot1, unsigned int OldNxtot2, int ObsBin) { 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); + c->SigmaTildeMuIndep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin), + OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); fastNLOTools::ExtendSigmaTildeX( - c->SigmaTildeMuRDep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); + c->SigmaTildeMuRDep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin), + OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); fastNLOTools::ExtendSigmaTildeX( - c->SigmaTildeMuFDep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); + c->SigmaTildeMuFDep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin), + OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); fastNLOTools::ExtendSigmaTildeX( - c->SigmaTildeMuRRDep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); + c->SigmaTildeMuRRDep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin), + OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); fastNLOTools::ExtendSigmaTildeX( - c->SigmaTildeMuFFDep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); + c->SigmaTildeMuFFDep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin), + OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); fastNLOTools::ExtendSigmaTildeX( - c->SigmaTildeMuRFDep[ObsBin], c->GetNxtot1(ObsBin), c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); + c->SigmaTildeMuRFDep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin), + OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue); } diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc index 973b7fdfe793770fbc36e67554874ccfdb6ef3b7..633578ea509d65723a1733b3ce067f6cd1c69025 100644 --- a/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc +++ b/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc @@ -581,26 +581,33 @@ namespace fastNLOTools { //______________________________________________________________________________ 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; - } + std::vector<T>& SigmaTildeX, unsigned int OldDimSize1, unsigned int NewDimSize1, + unsigned int OldDimSize2, unsigned int NewDimSize2, int NPDFDim, T InsertValue) { 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), InsertValue); - } + for (unsigned int i = OldDimSize1; i < NewDimSize1; i++) { + for (unsigned int j = 0; j <= i; j++) { + SigmaTildeX.insert(SigmaTildeX.begin() + (j * (j + 1) / 2), InsertValue); } - i++; } + } else if (NPDFDim == 2) { + for (unsigned int i = 0; i < OldDimSize2; i++) { + for (unsigned int j = OldDimSize1; j < NewDimSize1; j++) { + SigmaTildeX.insert(SigmaTildeX.begin() + i * NewDimSize1, InsertValue); + } + } + for (unsigned int i = OldDimSize2; i < NewDimSize2; i++) { + for (unsigned int j = 0; j < NewDimSize1; j++) { + SigmaTildeX.insert(SigmaTildeX.begin(), InsertValue); + } + } } else { - error["ExtendSigmaTildeX"] << "Unsupported NPdfDim: " << NPDFDim << endl; + error["ExtendSigmaTildeX"] << "Unsupported NPDFDim for x node density: " << NPDFDim << endl; 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); + template void ExtendSigmaTildeX<fastNLO::v1d>( + fastNLO::v2d&, unsigned int, unsigned int, unsigned int, unsigned int, int, fastNLO::v1d); + template void ExtendSigmaTildeX<fastNLO::v3d>( + fastNLO::v4d&, unsigned int, unsigned int, unsigned int, unsigned int, int, fastNLO::v3d); } // end namespace fastNLO diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCreate.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCreate.h index c6dfe6908f19018ba64ccc489bab33d76ab6b63f..af9b7b86243256744591ec9ab07cce0a1e38d221 100644 --- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCreate.h +++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCreate.h @@ -269,7 +269,7 @@ protected: void FlushCache(); //!< Fill weights from cache into table bool WarmupNeeded() const; - void ExtendAllFlexSigmaTilde(fastNLOCoeffAddFlex* c, int ObsBin); + void ExtendAllFlexSigmaTilde(fastNLOCoeffAddFlex* c, unsigned int OldNxtot1, unsigned int OldNxtot2, int ObsBin); struct fnloStats { diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h index b5c7f6930c34bc054b0843398f9bea3002cd6ded..5ffab5b7dab66930d8f2ac3bd77e13daa199ebcf 100644 --- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h +++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h @@ -95,7 +95,8 @@ namespace fastNLOTools { bool SameTails(fastNLO::v2d vector1, fastNLO::v2d vector2, double rtol=0.0); void ExtendHead(fastNLO::v1d& vector1, fastNLO::v1d& vector2); template <typename T> void ExtendSigmaTildeX( - std::vector<T>& SigmaTildeX, unsigned int DimSize1, unsigned int DimSize2, int NPDFDim, T InsertValue); + std::vector<T>& SigmaTildeX, unsigned int OldDimSize1, unsigned int NewDimSize1, + unsigned int OldDimSize2, unsigned int NewDimSize2, int NPDFDim, T InsertValue); };