diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
index da8b5dadf03be4ecbbfc772e0b3d139165f80220..ee1cbb30cf42beab9d79b76487a0cb671451823b 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
@@ -2446,116 +2446,8 @@ void fastNLOCreate::FillAllSubprocesses(const vector<vector<fnloEvent> >& events
//! events is expected to be of the form:
//! events[nscalevar][nsubproc]
- // KR TODO Check to replace wgtfac by _wo
- // cout << "WWW: wgtfac = " << wgtfac << "_wo = " << scen._wo << endl;
- const bool bFasterCode = false; // experimental developement: try to make code faster
- if (bFasterCode && !fIsWarmup && !fIsFlexibleScale) {
- // make filling code a little bit faster ... ~40%
- // if filling step "+=" is commented, then code is a factor ~6 faster
- fEvent = events[0][0];
- fScenario = scen;
-
- // KR: Also check for nan or inf in the faster code!!!
- if (!CheckWeightIsFinite()) return;
-
- const int ObsBin = GetBin();
- if (ObsBin < 0 || ObsBin >= (int)GetNObsBin()) return;
-
- // stats
- fStats._nEvPS++;
-
- fastNLOCoeffAddFix* c = (fastNLOCoeffAddFix*)GetTheCoeffTable();
- int NSubProc = c->GetNSubproc();
- int NPDFDim = c->GetNPDFDim();
- int NScaleVar = c->GetNScalevar();
- int NScaleNode = c->GetNScaleNode();
- // do interpolation
- 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;
- vector<pair<int,double> > nxup;
- unsigned int oldNxtot1 = c->XNode1[ObsBin].size();
- unsigned int oldNxtot2 = 0, newNxtot2 = 0;
- if (NPDFDim == 0) {
- if (xmin < xmax) {
- nxlo = fKernX1[ObsBin]->GetNodeValues(xmin);
- nxup = fKernX1[ObsBin]->GetNodeValues(xmax);
- } else {
- nxup = fKernX1[ObsBin]->GetNodeValues(xmax);
- nxlo = fKernX1[ObsBin]->GetNodeValues(xmin);
- }
- } else if (NPDFDim == 1) {
- nxlo = fKernX1[ObsBin]->GetNodeValues(xmin);
- nxup = fKernX1[ObsBin]->GetNodeValues(xmax);
- } else {
- nxlo = fKernX1[ObsBin]->GetNodeValues(xmin);
- nxup = fKernX2[ObsBin]->GetNodeValues(xmax);
- oldNxtot2 = c->XNode2[ObsBin].size();
- fastNLOTools::ExtendHead(c->XNode2[ObsBin], fKernX2[ObsBin]->fgrid);
- newNxtot2 = c->XNode2[ObsBin].size();
- }
- 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);
- fKernX2[ObsBin]->CheckX(xmax);
- ApplyPDFWeight(nxlo,xmin,fKernX1[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 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];
- for (unsigned int p = 0 ; p<events[is].size() ; p++) {
- double wgt = wgtfac * events[is][p]._w * nmu[m1].second / BinSize[ObsBin];
- // .......................................................................................
- for (unsigned int x1 = 0 ; x1<nxlo.size() ; x1++) {
- for (unsigned int x2 = 0 ; x2<nxup.size() ; x2++) {
- //int p = fEvent._p;
- int xminbin = nxlo[x1].first;
- int xmaxbin = nxup[x2].first;
- int proc = events[is][p]._p;
- HalfMatrixCheck(fEvent._x1,fEvent._x2,xminbin,xmaxbin,proc);
- int ixHM = GetXIndex(ObsBin,xminbin,xmaxbin);
- stm1[ixHM][proc] += nxlo[x1].second * nxup[x2].second * wgt;//nmu[m1].second * wp[p];
- // for (unsigned int x1 = 0 ; x1<nxup.size() ; x1++) {
- // for (unsigned int x2 = 0 ; x2<nxlo.size() ; x2++) {
- // int xmaxbin = nxup[x1].first;
- // int xminbin = nxlo[x2].first;
- // int proc = events[is][p]._p;
- // HalfMatrixCheck(fEvent._x1,fEvent._x2,xminbin,xmaxbin,proc);
- // int ixHM = GetXIndex(ObsBin,xminbin,xmaxbin);
- // // double ww = nxup[x1].second * nxlo[x2].second * wgt + stm1[ixHM][proc];
- // // double& sti = stm1[ixHM][proc];
- // // sti = ww; // this assignment is time consuming ?!
- // stm1[ixHM][proc] += nxup[x1].second * nxlo[x2].second * wgt;//nmu[m1].second * wp[p];
- // // c->SigmaTilde[ObsBin][is][nmu[m1].first][ixHM][proc] += nxup[x1].second * nxlo[x2].second * wgt;//nmu[m1].second * wp[p];
- }
- }
- }
- }
- }
- } else {
- for (unsigned int is = 0 ; is<events.size() ; is++) { // all scalevars
- FillAllSubprocesses(events[is], scen, is, wgtfac);
- }
+ for (unsigned int is = 0 ; is<events.size() ; is++) { // all scalevars
+ FillAllSubprocesses(events[is], scen, is, wgtfac);
}
}
@@ -2908,33 +2800,8 @@ void fastNLOCreate::FillContribution(int scalevar, const double wgtfac) {
exit(1);
}
- // hier
- // cout<<"Fill ObsBin="<<ObsBin<<", w="<<fEvent._w<<", x1="<<fEvent._x1<<", x2="<<fEvent._x2
- // <<", o0="<<fScenario._o[0]
- // <<", m1="<<fScenario._m1
- // <<", m2="<<fScenario._m2
- // <<", OB="<<fScenario._iOB
- // <<"\tproc="<<fEvent._p
- // <<", wr="<<fEvent._wr
- // <<", wf="<<fEvent._wf
- // <<endl;
- // static double wsum = 0;
- // wsum+= fEvent._w; //hier
- // cout<<" * wSum = "<<wsum<<endl;
-
-
- // --- statistics and weights**2
-
-
- //logger.debug["FillContributionFixHHC"]<<"The process Id is p = "<<p<<endl;
-
- //! read informatio from 'Event' and 'Scenario'
- //! do the interpolation
- //! and fill into the tables.
- //logger.debug["FillContributionFlexHHC"]<<endl;
-
if ( wgtfac != 1.0 ) {
- logger.error["FillContributionFlexHHC"]<<"Additional weight factor wgtfac = " << wgtfac << "not yet implemented for flex-scale tables, aborted!"<<endl;
+ logger.error[__func__]<<"Additional weight factor wgtfac = " << wgtfac << "not yet implemented for flex-scale tables, aborted!"<<endl;
exit(1);
}
@@ -3001,367 +2868,6 @@ void fastNLOCreate::FillContribution(int scalevar, const double wgtfac) {
//cout<<" * wSumW = "<<wsum<<endl;
}
-
-
-// ___________________________________________________________________________________________________
-void fastNLOCreate::FillContributionFixHHC(fastNLOCoeffAddFix* c, int ObsBin, int scalevar, const double wgtfac) {
- //! read informatio from 'Event' and 'Scenario'
- //! do the interpolation
- //! and fill into the tables.
- //logger.debug["FillContributionFixHHC"]<<endl;
-
- if ( wgtfac != 1.0 ) {
- logger.debug["FillContributionFixHHC"]<<"Experimental: Additional weight factor wgtfac = " << wgtfac << "not much tested so far!"<<endl;
- }
-
- if (fEvent._w == 0) return; // nothing todo.
-
- // do interpolation
- 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);
-
- const double mu = fScenario._m1 * fScaleFac[scalevar];
- const vector<pair<int,double> >& nmu = fKernMuS[ObsBin][scalevar]->GetNodeValues(mu);
-
- if (fApplyPDFReweight) {
- fKernX1[ObsBin]->CheckX(xmin);
- fKernX2[ObsBin]->CheckX(xmax);
- ApplyPDFWeight(nxlo,xmin,fKernX1[ObsBin]->GetGridPtr());
- ApplyPDFWeight(nxup,xmax,fKernX2[ObsBin]->GetGridPtr());
- }
-
-
- // fill grid
- if (!CheckWeightIsFinite()) return;
- double wgt = wgtfac * fEvent._w / BinSize[ObsBin];
- for (unsigned int x1 = 0 ; x1<nxlo.size() ; x1++) {
- for (unsigned int x2 = 0 ; x2<nxup.size() ; x2++) {
- int p = fEvent._p;
- int xminbin = nxlo[x1].first;
- int xmaxbin = nxup[x2].first;
- HalfMatrixCheck(fEvent._x1,fEvent._x2,xminbin,xmaxbin,p);
- //HalfMatrixCheck(xminbin,xmaxbin,p);
- int ixHM = GetXIndex(ObsBin,xminbin,xmaxbin);
- double w = wgt * nxlo[x1].second * nxup[x2].second;
- c->Fill(fEvent, ObsBin, ixHM, scalevar, nmu, nmu, p, w);
- }
- }
-}
-
-
-
-// ___________________________________________________________________________________________________
-void fastNLOCreate::FillContributionFlexHHC(fastNLOCoeffAddFlex* c, int ObsBin, const double wgtfac) {
- //! read informatio from 'Event' and 'Scenario'
- //! do the interpolation
- //! and fill into the tables.
- //logger.debug["FillContributionFlexHHC"]<<endl;
-
- if ( wgtfac != 1.0 ) {
- logger.error["FillContributionFlexHHC"]<<"Additional weight factor wgtfac = " << wgtfac << "not yet implemented for flex-scale tables, aborted!"<<endl;
- exit(1);
- }
-
- if (fEvent._w == 0 && fEvent._wf==0 && fEvent._wr==0 && fEvent._wrr==0 && fEvent._wff==0 && fEvent._wrf==0) return; // nothing todo.
- int NPDFDim = c->GetNPDFDim();
-
- // do interpolation
- //cout<<"try to interpol. ObsBin="<<ObsBin<<" ,x1="<<fEvent._x1<<", x2="<<fEvent._x2<<", mu1="<<Scenario._m1<<", mu2="<<Scenario._m2<<endl;
- 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> > 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);
- fastNLOTools::ExtendHead(c->XNode2[ObsBin], fKernX2[ObsBin]->fgrid);
- } else {
- nxup = fKernX1[ObsBin]->GetNodeValues(xmax);
- fastNLOTools::ExtendHead(c->XNode1[ObsBin], fKernX1[ObsBin]->fgrid);
- }
- ExtendAllFlexSigmaTilde(c, oldNxtot1, oldNxtot2, ObsBin);
-
- const vector<pair<int,double> >& nmu1 = fKernMu1[ObsBin]->GetNodeValues(fScenario._m1);
- const vector<pair<int,double> >& nmu2 = fKernMu2[ObsBin]->GetNodeValues(fScenario._m2);
-
- if (fApplyPDFReweight) {
- fKernX1[ObsBin]->CheckX(xmin);
- fKernX2[ObsBin]->CheckX(xmax);
- ApplyPDFWeight(nxlo,xmin,fKernX1[ObsBin]->GetGridPtr());
- if (c->NPDFDim > 1) {
- ApplyPDFWeight(nxup,xmax,fKernX2[ObsBin]->GetGridPtr());
- } else {
- ApplyPDFWeight(nxup,xmax,fKernX1[ObsBin]->GetGridPtr());
- }
- }
-
-
- // fill grid
- //if (!CheckWeightIsFinite()) return;
- for (unsigned int x1 = 0 ; x1<nxlo.size() ; x1++) {
- for (unsigned int x2 = 0 ; x2<nxup.size() ; x2++) {
- int xminbin = nxlo[x1].first;
- int xmaxbin = nxup[x2].first;
- int p = fEvent._p;
- HalfMatrixCheck(fEvent._x1,fEvent._x2,xminbin,xmaxbin,p);
- //HalfMatrixCheck(xminbin,xmaxbin,p);
- int ixHM = GetXIndex(ObsBin,xminbin,xmaxbin);
- double wfnlo = nxlo[x1].second * nxup[x2].second / BinSize[ObsBin];
- c->Fill(fEvent, ObsBin, ixHM, 0, nmu1, nmu2, p, wfnlo);
- }
- }
- //cout<<" * wSumW = "<<wsum<<endl;
-}
-
-
-
-// ___________________________________________________________________________________________________
-void fastNLOCreate::FillContributionFlexDIS(fastNLOCoeffAddFlex* c, int ObsBin, const double wgtfac) {
- //! read information from 'Event' and 'Scenario'
- //! do the interpolation
- //! and fill into the tables.
- //logger.debug["FillContributionFlexDIS"]<<endl;
-
- if ( wgtfac != 1.0 ) {
- logger.error["FillContributionFlexDIS"]<<"Additional weight factor wgtfac = " << wgtfac << "not yet implemented for flex-scale tables, aborted!"<<endl;
- exit(1);
- }
-
- if (fEvent._w == 0 && fEvent._wf==0 && fEvent._wr==0 && fEvent._wrr==0 && fEvent._wrf==0) return; // nothing todo.
-
- // do interpolation
- //cout<<"try to interpol. ObsBin="<<ObsBin<<" ,x1="<<fEvent._x1<<", x2="<<fEvent._x2<<", mu1="<<Scenario._m1<<", mu2="<<Scenario._m2<<endl;
-
- // todo, just: 'x'
- double x = fEvent._x1;
- vector<pair<int,double> > nx = fKernX1[ObsBin]->GetNodeValues(x);
- vector<pair<int,double> > nmu1 = fKernMu1[ObsBin]->GetNodeValues(fScenario._m1);
- vector<pair<int,double> > nmu2 = fKernMu2[ObsBin]->GetNodeValues(fScenario._m2);
-
-
- if (fApplyPDFReweight) {
- fKernX1[ObsBin]->CheckX(x);
- ApplyPDFWeight(nx,x,fKernX1[ObsBin]->GetGridPtr());
- }
-
- // fill grid
- if (!CheckWeightIsFinite()) return;
- /*
- for (unsigned int ix = 0 ; ix<nx.size() ; ix++) {
- int p = fEvent._p;
- int xIdx = nx[ix].first;
- //HalfMatrixCheck(xminbin,xmaxbin,p);
- //int ixHM = GetXIndex(ObsBin,xminbin,xmaxbin);
-
- for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) {
- for (unsigned int mu2 = 0 ; mu2<nmu2.size() ; mu2++) {
- double wfnlo = nx[ix].second * nmu1[m1].second * nmu2[mu2].second / BinSize[ObsBin];
- // if (! std::isfinite(wfnlo)) {
- // logger.error["FillContributionFlexDIS"]<<"Weight wfnlo is not finite, wfnlo = " << wfnlo << "!"<<endl;
- // logger.error["FillContributionFlexDIS"]<<"This should have been captured before, aborting ..."<<endl;
- // fKernX1[ObsBin]->PrintGrid();
- // fKernMu1[ObsBin]->PrintGrid();
- // fKernMu2[ObsBin]->PrintGrid();
- // cout<<"ix1="<<ix<<", im1="<<m1<<", im2="<<mu2<<endl;
- // cout<<"x1="<<nx[ix].second<<", ix="<<ix<<", xval="<<x<<endl;
- // cout<<"m1="<< nmu1[m1].second<<", m1="<<m1<<", mu1val="<<fScenario._m1<<endl;
- // cout<<"m2="<<nmu2[mu2].second<<", m2="<<mu2<<", mu2val="<<fScenario._m2<<endl;
- // exit(1);
- // }
- if (fEvent._w != 0) {
- //cout<<" Fill * : ix="<<ix<<", im1="<<nmu1[m1].first<<", im2="<<nmu2[mu2].first<<", p="<<p<<", w="<<fEvent._w * wfnlo<<endl;
- c->SigmaTildeMuIndep[ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._w * wfnlo;
- }
- if (fEvent._wf != 0) {
- //cout<<" Fill F : ix="<<ix<<", im1="<<nmu1[m1].first<<", im2="<<nmu2[mu2].first<<", p="<<p<<", w="<<fEvent._wf * wfnlo<<endl;
- c->SigmaTildeMuFDep [ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._wf * wfnlo;
- }
- if (fEvent._wr != 0) {
- //cout<<" Fill R : ix="<<ix<<", im1="<<nmu1[m1].first<<", im2="<<nmu2[mu2].first<<", p="<<p<<", w="<<fEvent._wr * wfnlo<<endl;
- c->SigmaTildeMuRDep [ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._wr * wfnlo;
- }
- if (fEvent._wrr != 0) {
- c->SigmaTildeMuRRDep [ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._wrr * wfnlo;
- }
- if (fEvent._wff != 0) {
- c->SigmaTildeMuFFDep [ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._wff * wfnlo;
- }
- if (fEvent._wrf != 0) {
- c->SigmaTildeMuRFDep [ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._wrf * wfnlo;
- }
- }
- }
- }
- */
-
- int p = fEvent._p;
- if (fEvent._w != 0) {
- for (unsigned int ix = 0 ; ix<nx.size() ; ix++) {
- int xIdx = nx[ix].first;
- for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) {
- for (unsigned int mu2 = 0 ; mu2<nmu2.size() ; mu2++) {
- double wfnlo = nx[ix].second * nmu1[m1].second * nmu2[mu2].second / BinSize[ObsBin];
- c->SigmaTildeMuIndep[ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._w * wfnlo;
- }
- }
- }
- }
- if (fEvent._wf != 0) {
- for (unsigned int ix = 0 ; ix<nx.size() ; ix++) {
- int xIdx = nx[ix].first;
- for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) {
- for (unsigned int mu2 = 0 ; mu2<nmu2.size() ; mu2++) {
- double wfnlo = nx[ix].second * nmu1[m1].second * nmu2[mu2].second / BinSize[ObsBin];
- c->SigmaTildeMuFDep [ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._wf * wfnlo;
- }
- }
- }
- }
- if (fEvent._wr != 0) {
- for (unsigned int ix = 0 ; ix<nx.size() ; ix++) {
- int xIdx = nx[ix].first;
- for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) {
- for (unsigned int mu2 = 0 ; mu2<nmu2.size() ; mu2++) {
- double wfnlo = nx[ix].second * nmu1[m1].second * nmu2[mu2].second / BinSize[ObsBin];
- c->SigmaTildeMuRDep [ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._wr * wfnlo;
- }
- }
- }
- }
- if (fEvent._wrr != 0) {
- for (unsigned int ix = 0 ; ix<nx.size() ; ix++) {
- int xIdx = nx[ix].first;
- for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) {
- for (unsigned int mu2 = 0 ; mu2<nmu2.size() ; mu2++) {
- double wfnlo = nx[ix].second * nmu1[m1].second * nmu2[mu2].second / BinSize[ObsBin];
- c->SigmaTildeMuRRDep [ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._wrr * wfnlo;
- }
- }
- }
- }
- if (fEvent._wff != 0) {
- for (unsigned int ix = 0 ; ix<nx.size() ; ix++) {
- int xIdx = nx[ix].first;
- for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) {
- for (unsigned int mu2 = 0 ; mu2<nmu2.size() ; mu2++) {
- double wfnlo = nx[ix].second * nmu1[m1].second * nmu2[mu2].second / BinSize[ObsBin];
- c->SigmaTildeMuFFDep [ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._wff * wfnlo;
- }
- }
- }
- }
- if (fEvent._wrf != 0) {
- for (unsigned int ix = 0 ; ix<nx.size() ; ix++) {
- int xIdx = nx[ix].first;
- for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) {
- for (unsigned int mu2 = 0 ; mu2<nmu2.size() ; mu2++) {
- double wfnlo = nx[ix].second * nmu1[m1].second * nmu2[mu2].second / BinSize[ObsBin];
- c->SigmaTildeMuRFDep [ObsBin][xIdx][nmu1[m1].first][nmu2[mu2].first][p] += fEvent._wrf * wfnlo;
- }
- }
- }
- }
-}
-
-
-// ___________________________________________________________________________________________________
-void fastNLOCreate::ExtendAllFlexSigmaTilde(
- fastNLOCoeffAddFlex* c, unsigned int OldNxtot1, unsigned int OldNxtot2, int ObsBin) {
- fastNLO::v3d insertValue = fastNLO::v3d(
- c->GetNScaleNode1(ObsBin),
- fastNLO::v2d(
- c->GetNScaleNode2(ObsBin),
- fastNLO::v1d(c->GetNSubproc())));
- unsigned int newNxtot1 = c->GetNxtot1(ObsBin);
- unsigned int newNxtot2 = c->GetNxtot2(ObsBin);
- fastNLOTools::ExtendSigmaTildeX(
- c->SigmaTildeMuIndep[ObsBin], OldNxtot1, newNxtot1,
- OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
- fastNLOTools::ExtendSigmaTildeX(
- c->SigmaTildeMuRDep[ObsBin], OldNxtot1, newNxtot1,
- OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
- fastNLOTools::ExtendSigmaTildeX(
- c->SigmaTildeMuFDep[ObsBin], OldNxtot1, newNxtot1,
- OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
- fastNLOTools::ExtendSigmaTildeX(
- c->SigmaTildeMuRRDep[ObsBin], OldNxtot1, newNxtot1,
- OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
- fastNLOTools::ExtendSigmaTildeX(
- c->SigmaTildeMuFFDep[ObsBin], OldNxtot1, newNxtot1,
- OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
- fastNLOTools::ExtendSigmaTildeX(
- c->SigmaTildeMuRFDep[ObsBin], OldNxtot1, newNxtot1,
- OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
-}
-
-
-// ___________________________________________________________________________________________________
-void fastNLOCreate::FillContributionFixDIS(fastNLOCoeffAddFix* c, int ObsBin, int scalevar, const double wgtfac) {
- //! read information from 'Event' and 'Scenario'
- //! do the interpolation
- //! and fill into the tables.
- //logger.debug["FillContributionFixDIS"]<<endl;
-
- if ( wgtfac != 1.0 ) {
- logger.warn["FillContributionFixDIS"]<<"Attention! Additional weight factor wgtfac = " << wgtfac << "not really tested so far!"<<endl;
- }
-
- if (fEvent._w == 0) return; // nothing todo.
- if (scalevar >= (int)fScaleFac.size()) {
- logger.error["FillContributionFixDIS"]<<"Error! Scale variation scalevar="<<scalevar<<" requested"
- <<", but only "<<fScaleFac.size()<<" variations are initialized. Exiting."<<endl;
- exit(3);
- }
-
- // do interpolation
- //cout<<"try to interpol. ObsBin="<<ObsBin<<" ,x1="<<fEvent._x1<<", x2="<<fEvent._x2<<", mu1="<<fScenario._m1<<endl;
-
- // todo, just: 'x'
- double x = fEvent._x1;
- vector<pair<int,double> > nx = fKernX1[ObsBin]->GetNodeValues(x);
- vector<pair<int,double> > nmu1 = fKernMuS[ObsBin][scalevar]->GetNodeValues(fScenario._m1);
-
-
- if (fApplyPDFReweight) {
- fKernX1[ObsBin]->CheckX(x);
- ApplyPDFWeight(nx,x,fKernX1[ObsBin]->GetGridPtr());
- }
-
- // fill grid
- if (!CheckWeightIsFinite()) return;
- for (unsigned int ix = 0 ; ix<nx.size() ; ix++) {
- int p = fEvent._p;
- int xIdx = nx[ix].first;
- //HalfMatrixCheck(xminbin,xmaxbin,p);
- //int ixHM = GetXIndex(ObsBin,xminbin,xmaxbin);
-
- for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) {
- double wfnlo = wgtfac * nx[ix].second * nmu1[m1].second / BinSize[ObsBin];
- if (! std::isfinite(wfnlo)) {
- logger.error["FillContributionFixDIS"]<<"Weight wfnlo is not finite, wfnlo = " << wfnlo << "!"<<endl;
- logger.error["FillContributionFixDIS"]<<"This should have been captured before, aborting ..."<<endl;
- fKernX1[ObsBin]->PrintGrid();
- fKernMu1[ObsBin]->PrintGrid();
- cout<<"ix1="<<ix<<", im1="<<m1<<endl;
- cout<<"x1="<<nx[ix].second<<", ix="<<ix<<", xval="<<x<<endl;
- cout<<"m1="<< nmu1[m1].second<<", m1="<<m1<<", mu1val="<<fScenario._m1<<endl;
- exit(1);
- }
- if (fEvent._w != 0) {
- //cout<<" Fill * : ix="<<xIdx<<", im1="<<nmu1[m1].first<<", im2="<<", p="<<p<<", w="<<fEvent._w * wfnlo<<endl;
- c->SigmaTilde[ObsBin][scalevar][nmu1[m1].first][xIdx][p] += fEvent._w * wfnlo;
- }
- }
- }
-}
-
-
-
// ___________________________________________________________________________________________________
void fastNLOCreate::FillRefContribution(int scalevar, const double wgtfac) {
//! This is a reference table.
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolBase.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolBase.cc
index 9a593a0bc6e7d2aea6ea59ccf4a04db4110aa5db..88fd0d71a3a494fadd01a6e8ba715f3d4e6f6be2 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolBase.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOInterpolBase.cc
@@ -34,28 +34,35 @@ fastNLOInterpolBase::fastNLOInterpolBase(double density, fastNLOGrid::GridType t
fExtendLow = true;
fvalmax = 1.0;
- double valMinH = -(nMinNodes - 1) / density;
+ double valMaxH;
+ const double valMinOffsetH = -(nMinNodes - 1) / density;
switch (fdm) {
case fastNLOGrid::kLinear:
- fvalmin = valMinH;
+ fvalmin = fvalmax - valMinOffsetH;
break;
case fastNLOGrid::kLogLog025:
- fvalmin = Function_loglog025_inv(valMinH);
+ valMaxH = Function_loglog025(fvalmax);
+ fvalmin = Function_loglog025_inv(valMaxH - valMinOffsetH);
break;
case fastNLOGrid::kLog10:
- fvalmin = Function_log10_inv(valMinH);
+ valMaxH = Function_log10(fvalmax);
+ fvalmin = Function_log10_inv(valMaxH - valMinOffsetH);
break;
case fastNLOGrid::kSqrtLog10:
- fvalmin = Function_sqrtlog10_inv(valMinH);
+ valMaxH = Function_sqrtlog10(fvalmax);
+ fvalmin = Function_sqrtlog10_inv(valMaxH - valMinOffsetH);
break;
case fastNLOGrid::kLogLog:
- fvalmin = Function_loglog_inv(valMinH);
+ valMaxH = Function_loglog(fvalmax);
+ fvalmin = Function_loglog_inv(valMaxH - valMinOffsetH);
break;
case fastNLOGrid::k3rdrtLog10:
- fvalmin = Function_3rdrtlog10_inv(valMinH);
+ valMaxH = Function_3rdrtlog10(fvalmax);
+ fvalmin = Function_3rdrtlog10_inv(valMaxH - valMinOffsetH);
break;
case fastNLOGrid::k4thrtLog10:
- fvalmin = Function_4thrtlog10_inv(valMinH);
+ valMaxH = Function_4thrtlog10(fvalmax);
+ fvalmin = Function_4thrtlog10_inv(valMaxH - valMinOffsetH);
break;
default:
error["MakeGridFromHGrid"]<<"Unknown grid type."<<endl;
@@ -384,7 +391,7 @@ bool fastNLOInterpolBase::CheckX(double& x) {
return true;
}
//printf("x=%e, %e; fgrid=%e, %e; ratio: %e\n",x,x-1.,fgrid[0],fgrid[0]-1,x/fgrid[0]-1);
- if ( x < fgrid[0] && false) {
+ if ( x < fgrid[0] && !fExtendLow) {
if ( x!=fLastVal[1] && fgrid[0]/x-1>1.e-6)
warn["CheckX"]<<"Value "<<x<<" is smaller than smallest node (min="<<fgrid[0]<<"). Using this first node."<<endl;
fLastVal[1] = x; // use this to monitor whenever there was an incident
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCreate.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCreate.h
index af9b7b86243256744591ec9ab07cce0a1e38d221..8bace7e196db4047d184f03fb82a362ac91aeabb 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCreate.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCreate.h
@@ -165,14 +165,6 @@ protected:
inline double CalcPDFReweight(double x) const;
//!< Fill contribution into table
void FillContribution(int scalevar = 0, const double wgtfac=1.0);
- //!< Fill flexible scale contribution in pp/ppbar
- void FillContributionFlexHHC(fastNLOCoeffAddFlex* c, int ObsBin, const double wgtfac=1.0);
- //!< fill flexible scale contribution in DIS
- void FillContributionFlexDIS(fastNLOCoeffAddFlex* c, int ObsBin, const double wgtfac=1.0);
- //!< fill fixed scale table in pp/ppbar
- void FillContributionFixHHC(fastNLOCoeffAddFix* c, int ObsBin, int scalevar, const double wgtfac=1.0);
- //!< fill fixed scale contribution in DIS
- void FillContributionFixDIS(fastNLOCoeffAddFix* c, int ObsBin, int scalevar, const double wgtfac=1.0);
//!< fill contribution if this is a reference table
void FillRefContribution(int scalevar = 0, const double wgtfac=1.0);
// KR: Don't see any use for shouldReadSteeringFile ==> remove it
@@ -269,8 +261,6 @@ protected:
void FlushCache(); //!< Fill weights from cache into table
bool WarmupNeeded() const;
- void ExtendAllFlexSigmaTilde(fastNLOCoeffAddFlex* c, unsigned int OldNxtot1, unsigned int OldNxtot2, int ObsBin);
-
struct fnloStats {
//! structre to keep track of statisics. Just for fun and information.