diff --git a/doc/tableformat/fastNLOTableFormatv25.ods b/doc/tableformat/fastNLOTableFormatv25.ods index 6cc3fa9ed593b2ccf1a8855c9e500902e665ab2c..53b4351ae83d36f825dd6ca4a7200f4e991cbe79 100644 Binary files a/doc/tableformat/fastNLOTableFormatv25.ods and b/doc/tableformat/fastNLOTableFormatv25.ods differ diff --git a/v2.3/generators/nlojet++/configure.ac b/v2.3/generators/nlojet++/configure.ac index 4a8472b07b3a572db499d77bd99107c5f1671759..0fcd5e524a0d47f94f99cd0da7b0368e85137d22 100644 --- a/v2.3/generators/nlojet++/configure.ac +++ b/v2.3/generators/nlojet++/configure.ac @@ -11,7 +11,7 @@ # Require minimal autoconf version of 2.69 is from 2012 AC_PREREQ([2.69]) # Define subproject fastNLO_interface_nlojet -AC_INIT([fastNLO_interface_nlojet], [2.3.1pre], [daniel.britzger@desy.de,klaus.rabbertz@cern.ch,g.sieber@cern.ch,stober@cern.ch,wobisch@fnal.gov]) +AC_INIT([fastNLO_interface_nlojet], [2.5.0], [daniel.britzger@desy.de,klaus.rabbertz@cern.ch,g.sieber@cern.ch,stober@cern.ch,wobisch@fnal.gov]) #Suppress verbose output when compiling (use make V=99 for verbose output) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) # Properly include subprojects diff --git a/v2.3/toolkit/data/check/InclusiveNJets_fnr0001midpHT_I723509_v23_fix-hhc-born-2jet_stat.log b/v2.3/toolkit/data/check/InclusiveNJets_fnr0001midpHT_I723509_v23_fix-hhc-born-2jet_stat.log index 99bdbd43a880e8aec55ba128953532f7dbc4ae5c..2e375907f8aad60d31c04c8d684129b57d1552fa 100644 --- a/v2.3/toolkit/data/check/InclusiveNJets_fnr0001midpHT_I723509_v23_fix-hhc-born-2jet_stat.log +++ b/v2.3/toolkit/data/check/InclusiveNJets_fnr0001midpHT_I723509_v23_fix-hhc-born-2jet_stat.log @@ -23,7 +23,7 @@ # Relative Statistical Uncertainties #--------------------------------------------------------------------------------- # bin cross section lower uncertainty upper uncertainty - #--------------------------------------------------------------------------------- + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1 8.01592840236E+05 -1.50893932747E-04 1.50893932747E-04 2 2.57832436144E+05 -1.19157110976E-04 1.19157110976E-04 3 6.90358759335E+04 -4.78717475807E-05 4.78717475807E-05 @@ -33,4 +33,4 @@ 7 1.84078315439E+02 -1.59591336380E-04 1.59591336380E-04 8 2.83497384445E+01 -8.90068444559E-05 8.90068444559E-05 9 3.04895884924E+00 -5.30297402284E-05 5.30297402284E-05 - #--------------------------------------------------------------------------------- + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/v2.3/toolkit/data/check/InclusiveNJets_fnr0001midpHT_I723509_v23_fix-hhc-nlo-2jet_stat.log b/v2.3/toolkit/data/check/InclusiveNJets_fnr0001midpHT_I723509_v23_fix-hhc-nlo-2jet_stat.log index a489c69c0e2739a9df9dfb1db76222185658532c..043d9f4b3ffa5e94e689513aaa5cd269b0a23d1d 100644 --- a/v2.3/toolkit/data/check/InclusiveNJets_fnr0001midpHT_I723509_v23_fix-hhc-nlo-2jet_stat.log +++ b/v2.3/toolkit/data/check/InclusiveNJets_fnr0001midpHT_I723509_v23_fix-hhc-nlo-2jet_stat.log @@ -23,7 +23,7 @@ # Relative Statistical Uncertainties #--------------------------------------------------------------------------------- # bin cross section lower uncertainty upper uncertainty - #--------------------------------------------------------------------------------- + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1 8.09746517361E+05 -8.17808102264E-03 8.17808102264E-03 2 2.51514173302E+05 -1.53654282706E-02 1.53654282706E-02 3 6.91227317135E+04 -1.05434381826E-03 1.05434381826E-03 @@ -33,4 +33,4 @@ 7 1.85065344888E+02 -6.31792468326E-03 6.31792468326E-03 8 2.82383223412E+01 -2.16054202482E-03 2.16054202482E-03 9 3.04174546319E+00 -2.14582701968E-03 2.14582701968E-03 - #--------------------------------------------------------------------------------- + #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/v2.3/toolkit/data/modify/SteerModify.str b/v2.3/toolkit/data/modify/SteerModify.str index a059e7ff3d36b9d11020d697ade98b3db083f1ef..973d309f505b7c30d58fae8d068046dbb87a1a94 100644 --- a/v2.3/toolkit/data/modify/SteerModify.str +++ b/v2.3/toolkit/data/modify/SteerModify.str @@ -108,14 +108,31 @@ ScDescript { # Example of final scenario description for CMS inclu #-------------------------------------------------------------------- # Add InfoBlocks with statistical/numerical uncertainty -# (For now only works with NNLOJET dat files!) #-------------------------------------------------------------------- -#InfoBlockFiles { # NNLOJET dat files for each order separately +#InfoBlockFlag1 0 # Default: 0 Statistical/numerical uncertainty +#InfoBlockFlag2 1 # Default: 1 = quadratic addition; else: 0 = linear addition + +# Either read uncertainty from files ... +#InfoBlockFiles { # For each order separately: NNLOJET dat files, fastNLO statunc log files, txt files using 1 or 2 columns # 2jet.LO.fnl3832_yb0_ys0_ptavgj12.dat # 2jet.NLO_only.fnl3832_yb0_ys0_ptavgj12.dat # 2jet.NNLO_only.fnl3832_yb0_ys0_ptavgj12.dat #} +#InfoBlockFileColumns { # Column(s) with rel. stat. uncertainty or with x section and abs. stat. uncertainty +# 4 +# 4 5 +#} + +# ... or provide them directly as columns inside this file +#InfoBlockStatUnc {{ +# dstrel_LO dstrel_NLO +# 0.001 0.01 +# 0.001 0.01 +# 0.001 0.01 +# 0.001 0.01 +#}} + #InfoBlockDescr { # Same description for each contribution # "Default statistical uncertainties from NNLOJET." #} diff --git a/v2.3/toolkit/fastnlotoolkit/fastNLOCoeffBase.cc b/v2.3/toolkit/fastnlotoolkit/fastNLOCoeffBase.cc index d7b54af53b84e1f709b6849b906f22fefd31bde0..97b90d68a5d6bb452168a8cb4847568815efd1b8 100644 --- a/v2.3/toolkit/fastnlotoolkit/fastNLOCoeffBase.cc +++ b/v2.3/toolkit/fastnlotoolkit/fastNLOCoeffBase.cc @@ -224,22 +224,39 @@ void fastNLOCoeffBase::MultiplyBin(unsigned int iObsIdx, double nfact) { // //________________________________________________________________________________________________________________ // Added to include CoeffInfoBlocks +// Check existence of InfoBlock of type, i.e. flags, (flag1,x); return error in case of multiple matches! bool fastNLOCoeffBase::HasCoeffInfoBlock(int fICoeffInfoBlockFlag1) const { bool result = false; for (int i=0; i<NCoeffInfoBlocks; i++) { - if ( ICoeffInfoBlockFlag1[i] == fICoeffInfoBlockFlag1 ) result = true; + if ( ICoeffInfoBlockFlag1[i] == fICoeffInfoBlockFlag1 ) { + if ( ! result ) { + result = true; + } else { + error["HasCoeffInfoBlocks"]<<"Aborted! Found multiple info blocks of type ICoeffInfoBlockFlag1 = "<<ICoeffInfoBlockFlag1[i]<<endl; + exit(201); + } + } } return result; } +// Check existence of InfoBlock of type, i.e. flags, (flag1,flag2); return error in case of multiple matches! bool fastNLOCoeffBase::HasCoeffInfoBlock(int fICoeffInfoBlockFlag1, int fICoeffInfoBlockFlag2) const { bool result = false; for (int i=0; i<NCoeffInfoBlocks; i++) { - if ( ICoeffInfoBlockFlag1[i] == fICoeffInfoBlockFlag1 && ICoeffInfoBlockFlag2[i] == fICoeffInfoBlockFlag2 ) result = true; + if ( ICoeffInfoBlockFlag1[i] == fICoeffInfoBlockFlag1 && ICoeffInfoBlockFlag2[i] == fICoeffInfoBlockFlag2 ) { + if ( ! result ) { + result = true; + } else { + error["HasCoeffInfoBlocks"]<<"Aborted! Found multiple info blocks of type ICoeffInfoBlockFlag1 = "<<ICoeffInfoBlockFlag1[i]<<endl; + exit(202); + } + } } return result; } +// Return first matching index; multiple blocks of same type, i.e. flags (flag1,x), are not allowed! int fastNLOCoeffBase::GetCoeffInfoBlockIndex(int fICoeffInfoBlockFlag1) { for (int i=0; i<NCoeffInfoBlocks; i++) { if ( ICoeffInfoBlockFlag1[i] == fICoeffInfoBlockFlag1 ) return i; @@ -247,6 +264,7 @@ int fastNLOCoeffBase::GetCoeffInfoBlockIndex(int fICoeffInfoBlockFlag1) { return -1; } +// Return first matching index; multiple blocks of same type, i.e. flags (flag1,flag2), are not allowed! int fastNLOCoeffBase::GetCoeffInfoBlockIndex(int fICoeffInfoBlockFlag1, int fICoeffInfoBlockFlag2) { for (int i=0; i<NCoeffInfoBlocks; i++) { if ( ICoeffInfoBlockFlag1[i] == fICoeffInfoBlockFlag1 && ICoeffInfoBlockFlag2[i] == fICoeffInfoBlockFlag2 ) return i; @@ -254,18 +272,25 @@ int fastNLOCoeffBase::GetCoeffInfoBlockIndex(int fICoeffInfoBlockFlag1, int fICo return -1; } -void fastNLOCoeffBase::AddCoeffInfoBlock(int fICoeffInfoBlockFlag1, int fICoeffInfoBlockFlag2, - std::vector<std::string> Description, std::string datfile) { +void fastNLOCoeffBase::AddCoeffInfoBlock(int fICoeffInfoBlockFlag1, int fICoeffInfoBlockFlag2, std::vector<std::string> Description, + std::vector<double> Uncertainty) { info["AddCoeffInfoBlocks"]<<"Adding additional InfoBlock with flags "<<fICoeffInfoBlockFlag1<<" and "<<fICoeffInfoBlockFlag2<<" to table contribution."<<endl; NCoeffInfoBlocks += 1; ICoeffInfoBlockFlag1.push_back(fICoeffInfoBlockFlag1); ICoeffInfoBlockFlag2.push_back(fICoeffInfoBlockFlag2); NCoeffInfoBlockDescr.push_back(Description.size()); CoeffInfoBlockDescript.push_back(Description); - std::vector<double> Uncertainty = fastNLOTools::ReadInfoBlockContent(datfile); + NCoeffInfoBlockCont.push_back(Uncertainty.size()); CoeffInfoBlockContent.push_back(Uncertainty); } +void fastNLOCoeffBase::AddCoeffInfoBlock(int fICoeffInfoBlockFlag1, int fICoeffInfoBlockFlag2, std::vector<std::string> Description, + std::string filename, unsigned int icola, unsigned int icolb) { + info["AddCoeffInfoBlocks"]<<"Adding additional InfoBlock reading data from file "<<filename<<endl; + std::vector<double> Uncertainty = fastNLOTools::ReadUncertaintyFromFile(filename, icola, icolb); + AddCoeffInfoBlock(fICoeffInfoBlockFlag1, fICoeffInfoBlockFlag2, Description, Uncertainty); +} + void fastNLOCoeffBase::ReadCoeffInfoBlocks(istream& table, int ITabVersionRead) { if (ITabVersionRead < 25000) { debug["ReadCoeffInfoBlocks"]<<"No additional info blocks allowed for table versions < 25000"<<endl; @@ -284,7 +309,7 @@ void fastNLOCoeffBase::ReadCoeffInfoBlocks(istream& table, int ITabVersionRead) } table >> iflag; ICoeffInfoBlockFlag2.push_back(iflag); - if (ICoeffInfoBlockFlag2[i] == 0) { + if (ICoeffInfoBlockFlag2[i] == 0 || ICoeffInfoBlockFlag2[i] == 1) { debug["ReadCoeffInfoBlocks"]<<"Found info block of type ICoeffInfoBlockFlag2 = "<<ICoeffInfoBlockFlag2[i]<<endl; } else { error["ReadCoeffInfoBlocks"]<<"Found info block of unknown type ICoeffInfoBlockFlag2 = "<<ICoeffInfoBlockFlag2[i]<<endl; diff --git a/v2.3/toolkit/fastnlotoolkit/fastNLODiffReader.cc b/v2.3/toolkit/fastnlotoolkit/fastNLODiffReader.cc index 0baa512bbab0a3a26c029eaf763327b9340f332b..209853851e634955c4f813acf9fd68878a786136 100644 --- a/v2.3/toolkit/fastnlotoolkit/fastNLODiffReader.cc +++ b/v2.3/toolkit/fastnlotoolkit/fastNLODiffReader.cc @@ -180,7 +180,7 @@ vector < double > fastNLODiffReader::GetDiffCrossSection() { // do the xpom integration logger.info["GetDiffCrossSection"]<<"Integrating xpom in "<<fxPoms.size()<<" slices. ["; fflush(stdout); - + fXSection_vs_xIPzIP.clear(); fXSection_vs_xIPzIP.resize(NObsBin); if ( fPrintxIPzIP ) { @@ -196,26 +196,26 @@ vector < double > fastNLODiffReader::GetDiffCrossSection() { if (i==0) logger.debug["GetDiffCrossSection"]<<"i="<<i<<"\tixp="<<ixp<<"\tfxpom="<<fxpom<<"\tXSection[i]="<<XSection[i]<<"\tfdxPoms[ixp]="<<fdxPoms[ixp]<<endl; xs[i] += XSection[i] * fdxPoms[ixp] ; xsLO[i] += XSection_LO[i] * fdxPoms[ixp] ; - - //vector<map<double,double> > fXSection_vsX1; - //std::vector < std::map< std::pair<double, double>, double > > fXSection_vs_xIPzIP; - if ( fPrintxIPzIP ) { - for ( auto xc : fXSection_vsX1[i] ) { - printf("%8d%14.6f%14.6f%14.6f\n",i,fxpom,xc.first,xc.second*fdxPoms[ixp]); - fXSection_vs_xIPzIP[i][make_pair(fxpom,xc.first)] = xc.second*fdxPoms[ixp]; - //cout<<i<<"\t"<<fxpom<<"\t"<<xc.first<<"\t"<<xc.second*fdxPoms[ixp]<<endl;; - } - } - // for ( auto xc : fXSection_vsX1[i] ) { - // fXSection_vs_xIPzIP[i][make_pair(xc.first,fxpom)] = xc.second*fdxPoms[ixp]; - // } - // for ( auto d : fXSection_vsX1[0] ) cout<<" | "<<d.first<<", "<<d.second; - // cout<<endl; + + //vector<map<double,double> > fXSection_vsX1; + //std::vector < std::map< std::pair<double, double>, double > > fXSection_vs_xIPzIP; + if ( fPrintxIPzIP ) { + for ( auto xc : fXSection_vsX1[i] ) { + printf("%8d%14.6f%14.6f%14.6f\n",i,fxpom,xc.first,xc.second*fdxPoms[ixp]); + fXSection_vs_xIPzIP[i][make_pair(fxpom,xc.first)] = xc.second*fdxPoms[ixp]; + //cout<<i<<"\t"<<fxpom<<"\t"<<xc.first<<"\t"<<xc.second*fdxPoms[ixp]<<endl;; + } + } + // for ( auto xc : fXSection_vsX1[i] ) { + // fXSection_vs_xIPzIP[i][make_pair(xc.first,fxpom)] = xc.second*fdxPoms[ixp]; + // } + // for ( auto d : fXSection_vsX1[0] ) cout<<" | "<<d.first<<", "<<d.second; + // cout<<endl; } // radek for(unsigned i = 0; i < NObsBin; ++i) { - //cout << "The basic cross section in bin "<<i <<" : " << XSection[i] << endl; + //cout << "The basic cross section in bin "<<i <<" : " << XSection[i] << endl; double sumMy=0; for(auto &x : fXSection_vsQ2[i]) { sumMy += x.second; @@ -230,16 +230,17 @@ vector < double > fastNLODiffReader::GetDiffCrossSection() { } logger.info>>"]"<<endl; logger.info["GetDiffCrossSection"]<< "Integrated interval in xpom: " << interv << endl; - + // set this cross section also to FastNLO mother class XSection = xs; XSection_LO = xsLO; - // k-factors - fastNLOReader::kFactor.resize(NObsBin); - for (unsigned int i = 0 ; i<NObsBin ; i++) { - fastNLOReader::kFactor[i] = fastNLOReader::XSection[i] / fastNLOReader::XSection_LO[i]; - } + // DEPRECATED + // // k-factors + // fastNLOReader::kFactor.resize(NObsBin); + // for (unsigned int i = 0 ; i<NObsBin ; i++) { + // fastNLOReader::kFactor[i] = fastNLOReader::XSection[i] / fastNLOReader::XSection_LO[i]; + // } return xs; } @@ -264,7 +265,7 @@ vector<double> fastNLODiffReader::GetXFX(double xp, double muf) const { int nx = -1; int nb = -1; for (int ib = 0 ; nb == -1 && ib<B_LO()->GetNObsBin() ; ib++) { - if (B_NLO() && (B_LO()->GetNxmax(ib) != B_NLO()->GetNxmax(ib))) + if (B_NLO() && (B_LO()->GetNxmax(ib) != B_NLO()->GetNxmax(ib))) logger.error["fastNLODiffReader::GetXFX"]<<"LO and NLO tables must have same number of x-bins."<<endl; for (int ix = 0 ; nx == -1 && ix<B_LO()->GetNxtot1(ib); ix++) { if (B_NLO() && (B_LO()->GetXNode1(ib,ix) != B_NLO()->GetXNode1(ib,ix))) diff --git a/v2.3/toolkit/fastnlotoolkit/fastNLOReader.cc b/v2.3/toolkit/fastnlotoolkit/fastNLOReader.cc index ba1e7aa6a3d922bf6ffdfe3b73d9cc6840880b01..aa090f458637951e84fa7a366ca0c50849ddb9b4 100644 --- a/v2.3/toolkit/fastnlotoolkit/fastNLOReader.cc +++ b/v2.3/toolkit/fastnlotoolkit/fastNLOReader.cc @@ -297,17 +297,6 @@ vector < double > xs = fnlo.GetCrossSection(); double* cs = &xs[0]; - Further you can access the "k-factor", which is calculated with all - 'contributions' that are switched on (e.g. non-perturbative corrections) - against the LO fixed-order contribution. - Remark: - - the proverbial k-factor is NLO vs. LO - - 1-loop threshold corrections are vs. LO - - 2-loop threshold corrections are vs. NLO - - non-perturbative corrections usually are vs. NLO - - vector < double > kFactors = fnlo.GetKFactors(); - 11. ---- Printing ---- // @@ -518,10 +507,6 @@ fastNLOReader::fastNLOReader(const fastNLOReader& other) : XSection(other.XSection), dXSection(other.dXSection), QScale(other.QScale), XSectionRef(other.XSectionRef), XSectionRefMixed(other.XSectionRefMixed), XSectionRef_s1(other.XSectionRef_s1), XSectionRef_s2(other.XSectionRef_s2) - // XSection_LO(other.XSection_LO), XSection(other.XSection), kFactor(other.kFactor), - // QScale_LO(other.QScale_LO), QScale(other.QScale), XSectionRef(other.XSectionRef), - // XSectionRefMixed(other.XSectionRefMixed), XSectionRef_s1(other.XSectionRef_s1), - // XSectionRef_s2(other.XSectionRef_s2) { //! Copy constructor OrderCoefficients(); // initialize pointers to fCoeff's @@ -908,7 +893,6 @@ bool fastNLOReader::SetContributionON(ESMCalculation eCalc , unsigned int Id , b } } - // if (!SetOld && SetOn && !fastNLOTools::IsEmptyVector(XSection_LO) ) { if (!SetOld && SetOn) { if (!c->GetIAddMultFlag()) { // if 'new' additive contribution, then refill PDF and alpha_s cache. // Fill alpha_s cache @@ -1026,13 +1010,6 @@ vector < vector < pair < int,int > > > fastNLOReader::GetSubprocIndices( const E } } -//______________________________________________________________________________ -vector < double > fastNLOReader::GetCrossSection() { - // Get fast calculated cross section - if (XSection.empty()) CalcCrossSection(); - return XSection; -} - //______________________________________________________________________________ vector < double > fastNLOReader::GetCrossSection(bool lNorm) { // Get fast calculated cross section @@ -1046,21 +1023,15 @@ vector < double > fastNLOReader::GetCrossSection(bool lNorm) { } //______________________________________________________________________________ -std::vector < std::map< double, double > > fastNLOReader::GetCrossSection_vs_x1() { - // Get fast calculated cross section - logger.warn<<"Function 'GetCrossSection_vs_x1' does _NOT_ return dSigma/dx but only the cross section contribution at the different x-nodes."<<endl; - logger.warn<<"In order to obtain dSigma/dx, the retured values must be divided by the step-size of the interpolation."<<endl; - if (XSection.empty()) CalcCrossSection(); - return fXSection_vsX1; -} - -//______________________________________________________________________________ -std::vector < std::map< double, double > > fastNLOReader::GetCrossSection_vs_x2() { - // Get fast calculated cross section - logger.warn<<"Function 'GetCrossSection_vs_x1' does _NOT_ return dSigma/dx but only the cross section contribution at the different x-nodes."<<endl; - logger.warn<<"In order to obtain dSigma/dx, the retured values must be divided by the step-size of the interpolation."<<endl; - if (XSection.empty()) CalcCrossSection(); - return fXSection_vsX2; +vector < double > fastNLOReader::GetUncertainty(bool lNorm) { + // Get uncertainty of fast calculated cross section stored in additional CoeffInfoBlocks + if (dXSection.empty()) CalcCrossSection(); + if (lNorm) { + logger.error["GetUncertainty"]<<"Additional uncertainty for normalised x sections not yet implemented; aborted!"<<endl; + exit(1); + } else { + return dXSection; + } } //______________________________________________________________________________ @@ -1133,6 +1104,24 @@ vector < double > fastNLOReader::GetNormCrossSection() { return XSectionNorm; } +//______________________________________________________________________________ +std::vector < std::map< double, double > > fastNLOReader::GetCrossSection_vs_x1() { + // Get fast calculated cross section + logger.warn<<"Function 'GetCrossSection_vs_x1' does _NOT_ return dSigma/dx but only the cross section contribution at the different x-nodes."<<endl; + logger.warn<<"In order to obtain dSigma/dx, the retured values must be divided by the step-size of the interpolation."<<endl; + if (XSection.empty()) CalcCrossSection(); + return fXSection_vsX1; +} + +//______________________________________________________________________________ +std::vector < std::map< double, double > > fastNLOReader::GetCrossSection_vs_x2() { + // Get fast calculated cross section + logger.warn<<"Function 'GetCrossSection_vs_x1' does _NOT_ return dSigma/dx but only the cross section contribution at the different x-nodes."<<endl; + logger.warn<<"In order to obtain dSigma/dx, the retured values must be divided by the step-size of the interpolation."<<endl; + if (XSection.empty()) CalcCrossSection(); + return fXSection_vsX2; +} + //______________________________________________________________________________ vector< vector < double > > fastNLOReader::GetCrossSection2Dim() { //! Get cross section as 2-dimensional vector according to defined binning @@ -1154,17 +1143,6 @@ vector< vector < double > > fastNLOReader::GetCrossSection2Dim() { } -//_DEPRECATED___________________________________________________________________ -vector < double > fastNLOReader::GetKFactors() { - // Get ratio of fast calculated NLO to LO cross section - // if (XSection.empty()) CalcCrossSection(); - logger.error["GetKFactors"]<<"This function is ambiguous and therefore deprecated, aborted!"<<endl; - logger.error["GetKFactors"]<<"Please derive K factors as ratios of the desired cross sections."<<endl; - exit(1); - return kFactor; -} - - //______________________________________________________________________________ vector < double > fastNLOReader::GetQScales() { // Get XSection weighted Q scale in bin @@ -1400,8 +1378,8 @@ void fastNLOReader::CalcCrossSection() { QScale[i] = QScale[i]/XSection[i]; } - // ---- Square root for statistical uncertainty combinaton ---- // - logger.debug["CalcCrossSection"]<<"Calculate stat. uncertainty from quadratic addition: sqrt(dXS)"<<endl; + // ---- Square root for summed statistical/numerical uncertainty ---- // + logger.debug["CalcCrossSection"]<<"Calculate statistical/numerical uncertainty from sqrt of summed contributions: sqrt(dXSection)"<<endl; for (unsigned int i=0; i<NObsBin; i++) { dXSection[i] = sqrt(dXSection[i]); } @@ -1503,9 +1481,9 @@ void fastNLOReader::CalcCrossSectionv21(fastNLOCoeffAddFlex* c) { if (!c) return; // Set up pointers to stored vectors - vector<double>* XS = &XSection; - vector<double>* dXS = &dXSection; - vector<double>* QS = &QScale; + vector<double>* XS = &XSection; + vector<double>* dXS = &dXSection; + vector<double>* QS = &QScale; // KR: Having different IXsectUnits in different contributions only works when // everything always scaled to Ipublunits (unique per table) @@ -1513,13 +1491,16 @@ void fastNLOReader::CalcCrossSectionv21(fastNLOCoeffAddFlex* c) { int xUnits = c->GetIXsectUnits(); logger.debug["CalcCrossSectionv21"]<<"Ipublunits = " << Ipublunits << ", xUnits = " << xUnits << endl; - // Check whether CoeffInfoBlock for relative statistical/numerical uncertainties (0,0) exists + // Check whether CoeffInfoBlock for statistical/numerical uncertainties (0,x) exists logger.debug["CalcCrossSectionv21"]<<"Checking on presence of statistical/numerical uncertainties ..."<<endl; - int iCIBIndex = c->GetCoeffInfoBlockIndex(0,0); + int iCIBIndex = -1; + int iCIBFlag2 = -1; std::vector < double > dCIBCont; - if ( iCIBIndex > -1 ) { + if ( c->HasCoeffInfoBlock(0) ) { + iCIBIndex = c->GetCoeffInfoBlockIndex(0); logger.debug["CalcCrossSectionv21"]<<"Found CoeffInfoBlock "<<iCIBIndex<<" with statistical/numerical uncertainties."<<endl; - dCIBCont = c->GetCoeffInfoContent(iCIBIndex); + iCIBFlag2 = c->GetCoeffInfoBlockFlag2(iCIBIndex); + dCIBCont = c->GetCoeffInfoContent(iCIBIndex); } else { logger.info["CalcCrossSectionv21"]<<"No CoeffInfoBlock found; uncertainties are initialised to zero."<<endl; } @@ -1575,10 +1556,16 @@ void fastNLOReader::CalcCrossSectionv21(fastNLOCoeffAddFlex* c) { if ( dCIBCont.empty() ) { dXS->at(i) += 0.; } else { - // Quadratical addition of independent statistical uncertainties - dXS->at(i) += XStmp*XStmp*dCIBCont[i]*dCIBCont[i]; - // Linear absolut addition of uncertainties (useful?) - // dXS->at(i) += fabs(XStmp)*dCIBCont[i]; + if ( iCIBFlag2 == 0 ) { + // Linear addition of absolute uncertainties; square is stored + dXS->at(i) = pow( (sqrt(dXS->at(i)) + fabs(XStmp)*dCIBCont[i]), 2 ); + } else if ( iCIBFlag2 == 1 ) { + // Quadratical addition of absolute uncertainties; square is stored + dXS->at(i) += XStmp*XStmp*dCIBCont[i]*dCIBCont[i]; + } else { + logger.error["CalcCrossSectionv21"]<<"Found illegal ICoeffInfoBlockFlag2 "<<iCIBFlag2<<", aborted!"<<endl; + exit(135); + } } } logger.debug["CalcCrossSectionv21"]<<"... leaving CalcCrossSectionv21."<<endl; @@ -1610,9 +1597,9 @@ void fastNLOReader::CalcCrossSectionv20(fastNLOCoeffAddFix* c) { } // Set up pointers to stored vectors - vector<double>* XS = &XSection; - vector<double>* dXS = &dXSection; - vector<double>* QS = &QScale; + vector<double>* XS = &XSection; + vector<double>* dXS = &dXSection; + vector<double>* QS = &QScale; // KR: Having different IXsectUnits in different contributions only works when // everything always scaled to Ipublunits (unique per table) @@ -1620,13 +1607,16 @@ void fastNLOReader::CalcCrossSectionv20(fastNLOCoeffAddFix* c) { int xUnits = c->GetIXsectUnits(); logger.debug["CalcCrossSectionv20"]<<"Ipublunits = " << Ipublunits << ", xUnits = " << xUnits << endl; - // Check whether CoeffInfoBlock for relative statistical/numerical uncertainties (0,0) exists + // Check whether CoeffInfoBlock for statistical/numerical uncertainties (0,x) exists logger.debug["CalcCrossSectionv20"]<<"Checking on presence of statistical/numerical uncertainties ..."<<endl; - int iCIBIndex = c->GetCoeffInfoBlockIndex(0,0); + int iCIBIndex = -1; + int iCIBFlag2 = -1; std::vector < double > dCIBCont; - if ( iCIBIndex > -1 ) { + if ( c->HasCoeffInfoBlock(0) ) { + iCIBIndex = c->GetCoeffInfoBlockIndex(0); logger.debug["CalcCrossSectionv20"]<<"Found CoeffInfoBlock "<<iCIBIndex<<" with statistical/numerical uncertainties."<<endl; - dCIBCont = c->GetCoeffInfoContent(iCIBIndex); + iCIBFlag2 = c->GetCoeffInfoBlockFlag2(iCIBIndex); + dCIBCont = c->GetCoeffInfoContent(iCIBIndex); } else { logger.info["CalcCrossSectionv20"]<<"No CoeffInfoBlock found; uncertainties are initialised to zero."<<endl; } @@ -1656,10 +1646,16 @@ void fastNLOReader::CalcCrossSectionv20(fastNLOCoeffAddFix* c) { if ( dCIBCont.empty() ) { dXS->at(i) += 0.; } else { - // Quadratical addition of independent statistical uncertainties - dXS->at(i) += XStmp*XStmp*dCIBCont[i]*dCIBCont[i]; - // Linear absolut addition of uncertainties (useful?) - // dXS->at(i) += fabs(XStmp)*dCIBCont[i]; + if ( iCIBFlag2 == 0 ) { + // Linear addition of absolute uncertainties + dXS->at(i) = pow( (sqrt(dXS->at(i)) + fabs(XStmp)*dCIBCont[i]), 2 ); + } else if ( iCIBFlag2 == 1 ) { + // Quadratical addition of absolute uncertainties + dXS->at(i) += XStmp*XStmp*dCIBCont[i]*dCIBCont[i]; + } else { + logger.error["CalcCrossSectionv20"]<<"Found illegal ICoeffInfoBlockFlag2 "<<iCIBFlag2<<", aborted!"<<endl; + exit(135); + } } } logger.debug["CalcCrossSectionv20"]<<"... leaving CalcCrossSectionv20."<<endl; @@ -1718,60 +1714,6 @@ void fastNLOReader::SetNewSqrtS(double newSqrtS, double SqrtStable) { } -//______________________________________________________________________________ -void fastNLOReader::SetCalculateSingleSubprocessOnly(int iSub) { - //! Calculate only a single subprocess with id=iSub - //! Please inform yourself before on the ordering of the subprocesses, - //! and if these are ordered identically for all contributions (LO,NLO,NNLO). - //! - //! iSub=-1 resets the calculation to all subprocesses - - /* - fSubprocActive.resize(13*13); - logger.info["SetCalculateSingleSubprocessOnly"]<<endl; - logger.info["SetCalculateSingleSubprocessOnly"]<<" *** Use this function carefully ***"<<endl; - logger.info["SetCalculateSingleSubprocessOnly"]<<" Please inform yourself about the meaning of the subprocess id for the given file."<<endl; - logger.info["SetCalculateSingleSubprocessOnly"]<<" This may also change for the different contribution (LO,NLO,NNLO) of a calculation."<<endl; - - if (iSub < 0) { - logger.info["SetCalculateSingleSubprocessOnly"]<<"Activating all contributions."<<endl; - fSubprocActive = std::vector<bool>(169,true); - } - else if ( iSub < 169 ) { - logger.info["SetCalculateSingleSubprocessOnly"]<<"Deactivating all contributions, but id="<<iSub<<endl; - fSubprocActive = std::vector<bool>(169,false); - fSubprocActive[iSub]=true; - } - */ - logger.warn["SetCalculateSingleSubprocessOnly"]<<endl<<"deprecated function, ignoring call"<<endl; -} - -//______________________________________________________________________________ -void fastNLOReader::SetCalculateSubprocesses( const std::vector<int>& iSub ){ - //! Calculate crosssection with several but not all Subprocesses included. - //! Information on the subprocesses available in the table can be retrieved - //! by calls to GetNSubproc and GetSubprocIndices. Also note, that subprocesses - //! cannot be split into single processes that were already merged at creation time - //! of the table. - //! - //! Use SetCalculateSingleSubprocessOnly( iSub=-1 ) to reset the calculation to use all processes. - - /* - fSubprocActive.resize(13*13); //this could probably dropped to the number returned byGetNSubproc - logger.info["SetCalculateSubprocesses"]<<endl; - logger.info["SetCalculateSubprocesses"]<<" *** Use this function carefully ***"<<endl; - logger.info["SetCalculateSubprocesses"]<<" Please inform yourself about the meaning of the subprocess id for the given file."<<endl; - logger.info["SetCalculateSubprocesses"]<<" This may also change for the different contribution (LO,NLO,NNLO) of a calculation."<<endl; - - fSubprocActive = std::vector<bool>(13*13,false); - for ( int i = 0; i < iSub.size(); i++ ) - if ( iSub[i] < 13*13 ) - fSubprocActive[iSub[i]] = true; - */ - logger.warn["SetCalculateSubprocesses"]<<endl<<"deprecated function, ignoring call"<<endl; - -} - //______________________________________________________________________________ void fastNLOReader::SelectProcesses( const std::vector< std::pair<int,int> >& proclist ) { //! Selects subprocesses given in proclist. proclist is a vector of pairs each identifying @@ -3293,20 +3235,6 @@ int fastNLOReader::ContrId(const ESMCalculation eCalc, const ESMOrder eOrder) co //______________________________________________________________________________ -//_DEPRECATED___________________________________________________________________ -void fastNLOReader::PrintTableInfo(const int iprint) const { - logger.warn["PrintTableInfo"]<<"This function is deprecated!"<<endl; - logger.warn["PrintTableInfo"]<<"Please use PrintContributionSummary instead."<<endl; -} - - -//_DEPRECATED___________________________________________________________________ -void fastNLOReader::PrintFastNLOTableConstants(const int iprint) const { - logger.warn["PrintTableInfo"]<<"This function is deprecated!"<<endl; - logger.warn["PrintTableInfo"]<<"Please use Print instead."<<endl; -} - - //______________________________________________________________________________ void fastNLOReader::PrintContributionSummary(int iprint) const { //! this function is inherited from fastNLOTable. @@ -3358,7 +3286,6 @@ void fastNLOReader::PrintCrossSections() const { printf(" # ----> from %9.3f to %9.3f in %s <----\n",GetObsBinLoBound(i,0),GetObsBinUpBound(i,0),GetDimLabel(0).c_str()); lobindim2 = GetObsBinLoBound(i,0); } - // printf(" # %4.0f | %9.3f - %9.3f % 9.4e % 5.2f |\n",i*1.,GetObsBinLoBound(i,1),GetObsBinUpBound(i,1),xs[i],kFactor[i]); printf(" # %4.0f | %9.3f - %9.3f % 9.4e |\n",i*1.,GetObsBinLoBound(i,1),GetObsBinUpBound(i,1),xs[i]); } } @@ -3441,8 +3368,6 @@ void fastNLOReader::PrintCrossSectionsWithReference() { printf(" # ----> from %9.3f to %9.3f in %s <----\n",GetObsBinLoBound(i,0),GetObsBinUpBound(i,0),GetDimLabel(0).c_str()); lobindim2 = GetObsBinLoBound(i,0); } - // printf(" # %4.0f | %9.3f - %9.3f % 9.4e % 5.3f | % 9.4e % 5.4f\n", - // i*1.,GetObsBinLoBound(i,1),GetObsBinUpBound(i,1),xs[i],kFactor[i],xsref[i],(xs[i]-xsref[i])/xsref[i]*100.); printf(" # %4.0f | %9.3f - %9.3f % 9.4e | % 9.4e % 5.4f\n", i*1.,GetObsBinLoBound(i,1),GetObsBinUpBound(i,1),xs[i],xsref[i],(xs[i]-xsref[i])/xsref[i]*100.); } @@ -3459,31 +3384,6 @@ void fastNLOReader::PrintCrossSectionsWithReference() { } -//_DEPRECATED___________________________________________________________________ -void fastNLOReader::PrintCrossSectionsData() const { - logger.warn["PrintCrossSectionsData"]<<"This function is deprecated!"<<endl; - logger.warn["PrintCrossSectionsData"]<<"Please check fnlo-tk-cppread.cc for default print out."<<endl; - return; -} - - -//_DEPRECATED___________________________________________________________________ -void fastNLOReader::PrintCrossSectionsDefault(const vector <double> kthc) const { - logger.warn["PrintCrossSectionsDefault"]<<"This function is deprecated!"<<endl; - logger.warn["PrintCrossSectionsDefault"]<<"Please check fnlo-tk-cppread.cc for default print out."<<endl; - exit(1); - return; -} - - -//_DEPRECATED___________________________________________________________________ -void fastNLOReader::RunFastNLODemo() { - logger.warn["RunFastNLODemo"]<<"This function is deprecated!"<<endl; - logger.warn["RunFastNLODemo"]<<"Please check fnlo-tk-examples.cc for examples."<<endl; - return; -} - - //______________________________________________________________________________ double fastNLOReader::RescaleCrossSectionUnits(double binsize, int xunits) { //! @@ -3503,13 +3403,6 @@ double fastNLOReader::RescaleCrossSectionUnits(double binsize, int xunits) { } -//______________________________________________________________________________ -XsUncertainty fastNLOReader::GetScaleUncertainty(const EScaleUncertaintyStyle eScaleUnc) { - XsUncertainty XsUnc = GetScaleUncertainty(eScaleUnc, false); - return XsUnc; -} - - //______________________________________________________________________________ XsUncertainty fastNLOReader::GetScaleUncertainty(const EScaleUncertaintyStyle eScaleUnc, bool lNorm) { // Get 2- or 6-point scale uncertainty @@ -3587,34 +3480,6 @@ std::vector< std::vector<double> > fastNLOReader::GetScaleUncertaintyVec(const E // Added to include CoeffInfoBlocks // -//______________________________________________________________________________ -vector < double > fastNLOReader::GetUncertainty() { - // Get uncertainty of fast calculated cross section stored in additional CoeffInfoBlocks - if (dXSection.empty()) CalcCrossSection(); - return dXSection; -} - -//______________________________________________________________________________ -vector < double > fastNLOReader::GetUncertainty(bool lNorm) { - // Get uncertainty of fast calculated cross section stored in additional CoeffInfoBlocks - if (dXSection.empty()) CalcCrossSection(); - if (lNorm) { - // vector < double > XNorm = GetNormCrossSection(); - // return XNorm; - logger.error["GetUncertainty"]<<"Additional uncertainty for normalised x sections not yet implemented; aborted!"<<endl; - exit(1); - } else { - return dXSection; - } -} - -//________________________________________________________________________________________________________________ -XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUnc) { - XsUncertainty XsUnc = GetAddUncertainty(eAddUnc, false); - return XsUnc; -} - - //______________________________________________________________________________ XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUnc, bool lNorm) { // @@ -3637,7 +3502,7 @@ XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUn XsUnc.dxsl.push_back(0); } } else if (eAddUnc == kAddStat) { - logger.info["GetAddUncertainty"]<<"Statistical uncertainties selected."<<endl; + logger.info["GetAddUncertainty"]<<"Statistical/numerical uncertainties selected."<<endl; for (unsigned int iobs = 0; iobs < NObsBin; iobs++) { XsUnc.xs.push_back(MyXSection[iobs]); XsUnc.dxsu.push_back(MydXSection[iobs]); @@ -3649,7 +3514,6 @@ XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUn exit(1); } - //! Divide by cross section != 0 to give relative uncertainties for (unsigned int iobs = 0; iobs < NObsBin; iobs++) { if (fabs(XsUnc.xs[iobs]) > DBL_MIN) { diff --git a/v2.3/toolkit/fastnlotoolkit/fastNLOTools.cc b/v2.3/toolkit/fastnlotoolkit/fastNLOTools.cc index fea7ffad1158258b53a90f6c944234200c56c9d5..7947823dd9a3f86e25733b2dd02fe8cbecfc69d0 100644 --- a/v2.3/toolkit/fastnlotoolkit/fastNLOTools.cc +++ b/v2.3/toolkit/fastnlotoolkit/fastNLOTools.cc @@ -392,26 +392,39 @@ namespace fastNLOTools { } //______________________________________________________________________________ - std::vector <double> ReadInfoBlockContent(std::string filename) { + std::vector <double> ReadUncertaintyFromFile(std::string filename, unsigned int icola, unsigned int icolb) { std::string extension = ""; std::ifstream infile; std::string line; std::vector <double> Uncertainty; - //! Determine extension to differentiate for parsing among - // - NLOJet++ --> fnlo-tk-statunc .log files - // - NNLOJET --> .dat files + //! Determine extension to differentiate for parsing + //! - fnlo-tk-statunc: 'log' file extension; column numbers not needed, rel. stat. uncertainty = col #4 + //! - NNLOJET dat file: 'dat' file extension; column numbers not needed, rel. stat. uncertainty = (col #5 / col #4) + //! - Generic txt file: 'txt' file extension; only icola --> rel. stat. uncertainty = col #icola + //! - icol a & b --> rel. stat. uncertainty = col #icolb / #icola if ( filename.find_last_of(".") != std::string::npos ) { extension = filename.substr(filename.find_last_of(".")+1); } - info["ReadInfoBlockContent"]<<"Reading additional InfoBlock content from file: " << filename <<endl; + if ( extension != "dat" && extension != "log" && extension != "txt" ) { + error["ReadUncertaintyFromFile"]<<"Unknown filename extension, aborted! filename = " << filename <<endl; + exit(34); + } else if ( extension == "txt" && icola == 0) { + error["ReadUncertaintyFromFile"]<<"'txt' file found, but column specification is missing, aborted! icola " << icola <<endl; + exit(35); + } else if ( extension == "txt" && (icola > 10 || icolb > 10) ) { + error["ReadUncertaintyFromFile"]<<"'txt' file found, but column specification is too large, aborted! icola, icolb = " << icola << ", " << icolb <<endl; + exit(35); + } else { + info["ReadUncertaintyFromFile"]<<"Reading additional uncertainty content from file: " << filename <<endl; + } + infile.open(filename); if (infile.is_open()) { int iline = 0; bool lline = false; // Read line-by-line while(std::getline(infile, line)) { - cout << "line = " << line << endl; // Put line into stringstream and read word-by-word std::istringstream iss(line); std::string word, word1, word2; @@ -443,24 +456,40 @@ namespace fastNLOTools { if ( word == "#-" ) { lline = ! lline; } else if ( lline ) { - cout << "AAAAA word = " << word << endl; // Skip second & third word of each uncertainty line (x section; lower uncertainty) iss >> word; iss >> word; double dxsrel; iss >> dxsrel; - cout << "dxsrel = " << dxsrel << endl; Uncertainty.push_back(dxsrel); iline += 1; - cout << "BBBBB iline = " << iline << endl; + } + } + // For 'txt' extension either read column #icola or divide column #icolb / #icola; max col = 10 + else if ( extension == "txt" ) { + double a = 0; + double b = 0; + for ( unsigned int ic = 1; ic<11; ic++ ) { + if ( ic == icola ) a = std::stod(word); + if ( ic == icolb ) b = std::stod(word); + iss >> word; + } + if ( icolb == 0 ) { + Uncertainty.push_back(a); + } else { + if ( fabs(a) > DBL_MIN ) { + Uncertainty.push_back(b/a); + } else { + Uncertainty.push_back(0); + } } } else { - error["ReadInfoBlockContent"]<<"Unkown file extension, aborted! Filename is: " << filename <<endl; - exit(32); + error["ReadUncertaintyFromFile"]<<"Unknown filename extension, aborted! filename = " << filename <<endl; + exit(34); } } } else { - error["ReadInfoBlockContent"]<<"Cannot read InfoBlock content, aborted! Filename is: " << filename <<endl; + error["ReadUncertaintyFromFile"]<<"Cannot read from file, aborted! filename is: " << filename <<endl; exit(33); } return Uncertainty; diff --git a/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffBase.h b/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffBase.h index 9ac3f2490f9996d00c36d57f005ac7b63e4eb764..678cd6b7bf4cf9c496a1500978d0e9135aac8529 100644 --- a/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffBase.h +++ b/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffBase.h @@ -85,9 +85,18 @@ public: bool HasCoeffInfoBlock(int ICoeffInfoBlockFlag1, int ICoeffInfoBlockFlag2) const; int GetCoeffInfoBlockIndex(int ICoeffInfoBlockFlag1); int GetCoeffInfoBlockIndex(int ICoeffInfoBlockFlag1, int ICoeffInfoBlockFlag2); + int GetCoeffInfoBlockFlag1(int Index) const { return ICoeffInfoBlockFlag1[Index]; }; + int GetCoeffInfoBlockFlag2(int Index) const { return ICoeffInfoBlockFlag2[Index]; }; + void SetCoeffInfoBlockFlag1(int Index, int iFlag1) { ICoeffInfoBlockFlag1[Index] = iFlag1; }; + void SetCoeffInfoBlockFlag2(int Index, int iFlag2) { ICoeffInfoBlockFlag2[Index] = iFlag2; }; std::vector < double > GetCoeffInfoContent(int Index) const { return CoeffInfoBlockContent[Index]; }; int GetNCoeffInfoBlocks() const {return NCoeffInfoBlocks;} - void AddCoeffInfoBlock(int ICoeffInfoBlockFlag1, int ICoeffInfoBlockFlag2, std::vector<std::string> Description, std::string datfile); + // Provide uncertainty via input vector + void AddCoeffInfoBlock(int ICoeffInfoBlockFlag1, int ICoeffInfoBlockFlag2, std::vector<std::string> Description, + std::vector<double> uncertainty); + // Provide uncertainty reading from filename + void AddCoeffInfoBlock(int ICoeffInfoBlockFlag1, int ICoeffInfoBlockFlag2, std::vector<std::string> Description, + std::string filename, unsigned int icola = 0, unsigned int icolb = 0); protected: void ReadBase(std::istream& table, int ITabVersionRead); diff --git a/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in b/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in index 3493edfbd1018ab956c93d34994db7e5e96d17f9..1ff6eeb7be2f8c252e69230e0b1bafe043216fb8 100644 --- a/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in +++ b/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in @@ -112,7 +112,7 @@ namespace fastNLO { enum EAddUncertaintyStyle { kAddNone = 0, // no additional uncertainty - kAddStat = 1 // statistical uncertainty ( i.e. uncorrelated and symmetric) + kAddStat = 1 // statistical/numerical uncertainty }; enum EPDFUncertaintyStyle { diff --git a/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h b/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h index 12e3d7c2dc4836fe44ff7abaedc242d5657a94cf..7a1e114673b821ba4cd47c99293d39a3a0e33f92 100644 --- a/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h +++ b/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h @@ -18,7 +18,7 @@ class fastNLOReader : public fastNLOTable , public fastNLOPDFLinearCombinations //! //! fastNLOReader. //! Abstract base class for evaluation of fastNLO tables. - //! Instantiatians must implement functions for PDF and alpha_s access. + //! Instantiations must implement functions for PDF and alpha_s access. //! public: @@ -51,8 +51,6 @@ public: void SelectProcesses( const std::string& processes, bool symmetric = true ); //!< tries to select the specified subprocesses for calculation. Prints a warning on failure. // ---- setters for specific options ---- // - void SetCalculateSingleSubprocessOnly(int iSub); //!< only use one subprocess for calculation - void SetCalculateSubprocesses( const std::vector< int >& iSub ); //!< use several but not all subprocesses in calculation void SetNewSqrtS(double NewSqrtS, double OldSqrtS=0 ); // ---- setters for scales of MuVar tables ---- // @@ -86,25 +84,20 @@ public: void CalcCrossSection(); double RescaleCrossSectionUnits(double binsize, int xunits); // Rescale according to kAbsoluteUnits and Ipublunits settings - std::vector < double > GetCrossSection(); //!< Return vector with all cross section values - std::vector < double > GetCrossSection(bool lNorm); //!< Return vector with all cross section values, normalize on request - std::vector < double > GetNormCrossSection(); //!< Return vector with all normalized cross section values - std::vector < double > GetUncertainty(); //!< Return vector with additional uncertainty of cross section values - std::vector < double > GetUncertainty(bool lNorm); //!< Return vector with additional uncertainty of cross section values, normalise on request (NOT YET IMPLEMENTED) + std::vector < double > GetCrossSection(bool lNorm = false); //!< Return vector with all cross section values, normalize on request + std::vector < double > GetUncertainty(bool lNorm = false); //!< Return vector with additional uncertainty of cross section values, normalise on request (NOT YET IMPLEMENTED) + std::vector < double > GetNormCrossSection(); //!< Return vector with all normalized cross section values std::vector < std::map< double, double > > GetCrossSection_vs_x1(); //! Cross section vs. x1 ( XSection_vsX1[bin][<x,xs>] ) std::vector < std::map< double, double > > GetCrossSection_vs_x2(); //! Cross section vs. x2 ( XSection_vsX1[bin][<x,xs>] ) std::vector < double > GetReferenceCrossSection(); - std::vector < double > GetKFactors(); ///< DEPRECATED std::vector < double > GetQScales(); //!< Order (power of alpha_s) rel. to LO: 0 --> LO, 1 --> NLO std::vector < std::vector < double > > GetCrossSection2Dim(); //! Return struct with vectors containing the cross section values and the selected uncertainty - XsUncertainty GetScaleUncertainty( const fastNLO::EScaleUncertaintyStyle eScaleUnc ); - XsUncertainty GetScaleUncertainty( const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm ); - XsUncertainty GetAddUncertainty( const fastNLO::EAddUncertaintyStyle eAddUnc ); - XsUncertainty GetAddUncertainty( const fastNLO::EAddUncertaintyStyle eAddUnc, bool lNorm ); + XsUncertainty GetScaleUncertainty( const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false); + XsUncertainty GetAddUncertainty( const fastNLO::EAddUncertaintyStyle eAddUnc, bool lNorm = false); //! Function for use with pyext (TODO: Clean this up) std::vector< std::vector<double> > GetScaleUncertaintyVec( const fastNLO::EScaleUncertaintyStyle eScaleUnc ); std::vector< std::vector<double> > GetAddUncertaintyVec( const fastNLO::EAddUncertaintyStyle eAddUnc ); @@ -136,13 +129,6 @@ public: void PrintCrossSections() const; //!< Print cross sections (optimized for double-differential tables) void PrintCrossSectionsWithReference(); - // DEPRECATED - void PrintTableInfo(const int iprint = 0) const; // DEPRECATED, use PrintContributionSummary instead - void PrintFastNLOTableConstants(const int iprint = 2) const; // DEPRECATED, use Print instead - void PrintCrossSectionsData() const; ///< DEPRECATED - void PrintCrossSectionsDefault(std::vector<double> kthc = std::vector<double>()) const; ///< DEPRECATED - void RunFastNLODemo(); ///< DEPRECATED - // ---- Test virtual functions for reasonable values. ---- // bool TestXFX(); //!< Test if XFX reasonable values bool TestAlphas(); //!< Test if EvolvaAlphas returns a reasonable value @@ -255,8 +241,8 @@ protected: // ---- Cross sections ---- // std::vector < double > XSection_LO; std::vector < double > XSection; - std::vector < double > dXSection; // stored uncertainty of x section - std::vector < double > kFactor; + std::vector < double > dXSection; // uncertainty sum of x section + std::vector < double > dX2Section; // squared uncertainty sum of x section std::vector < double > QScale_LO; std::vector < double > QScale; std::vector < std::map< double, double > > fXSection_vsX1; //! Cross section vs. x ( XSection_vsX1[bin][<x,xs>] ) diff --git a/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h b/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h index 39e528fd0117ba048e848709c3865f29a66c57dc..8cb5e0996e42a3456e5325cf130bcb50448c4119 100644 --- a/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h +++ b/v2.3/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h @@ -80,7 +80,12 @@ namespace fastNLOTools { bool ReadMagicNo(std::istream& table); //!< Read and check magic number from table. void PutBackMagicNo(std::istream& table); //!< Reset magic number, such that it can be recognized by other reading routines - std::vector <double> ReadInfoBlockContent(std::string filename); + //! Parse filename for uncertainties + //! - fnlo-tk-statunc: 'log' file extension; column numbers not needed, rel. stat. uncertainty = col #4 + //! - NNLOJET dat file: 'dat' file extension; column numbers not needed, rel. stat. uncertainty = (col #5 / col #4) + //! - Generic txt file: 'txt' file extension; only icola --> rel. stat. uncertainty = col #icola + //! - icol a & b --> rel. stat. uncertainty = col #icolb / #icola + std::vector <double> ReadUncertaintyFromFile(std::string filename, unsigned int icola = 0, unsigned int icolb = 0); }; diff --git a/v2.3/toolkit/src/fnlo-tk-modify.cc b/v2.3/toolkit/src/fnlo-tk-modify.cc index 19699965b45ba4de66c15f569f4cfaf0673e0577..d858dfab1480b84350c95243f77420d6c7caf97b 100644 --- a/v2.3/toolkit/src/fnlo-tk-modify.cc +++ b/v2.3/toolkit/src/fnlo-tk-modify.cc @@ -181,6 +181,7 @@ int main(int argc, char** argv) { //! --- Print input table information info["fnlo-tk-modify"] << "Trying to read input table " << STRING(InTable) << endl; fastNLOTable table(CHAR(InTable)); + unsigned int nobs = table.GetNObsBin(); if ( EXIST(PrintInputA1) ) { if ( BOOL(PrintInputA1) ) table.PrintHeader(1); } @@ -244,7 +245,6 @@ int main(int argc, char** argv) { if ( EXIST(BinSizeFactor) ) { double fac = DOUBLE(BinSizeFactor); info["fnlo-tk-modify"]<<"Multiplying all bin sizes by factor " << fac << "!" << endl; - unsigned int nobs = table.GetNObsBin(); for (unsigned int iObs=0; iObs<nobs; iObs++) { table.MultiplyBinSize(iObs,fac); } @@ -253,7 +253,6 @@ int main(int argc, char** argv) { if ( !DOUBLE_ARR(BinSize).empty() ) { vector<double> fac = DOUBLE_ARR(BinSize); info["fnlo-tk-modify"]<<"Multiplying bin sizes by provided factors!"<<endl; - unsigned int nobs = table.GetNObsBin(); if ( nobs != fac.size() ) { error["fnlo-tk-modify"]<<"You need the same number of multiplicative factors, nfact = " << fac.size() << ", than bins in the table, nobsbin = " << nobs << ". Aborted!" << endl; @@ -268,7 +267,6 @@ int main(int argc, char** argv) { if ( !DOUBLE_ARR(MultCoeff).empty() ) { vector<double> fac = DOUBLE_ARR(MultCoeff); info["fnlo-tk-modify"]<<"Multiplying by provided factors all coefficients of additive contributions to observable bins!"<<endl; - unsigned int nobs = table.GetNObsBin(); if ( nobs != fac.size() ) { error["fnlo-tk-modify"]<<"You need the same number of multiplicative factors, nfact = " << fac.size() << ", than bins in the table, nobsbin = " << nobs << ". Aborted!" << endl; @@ -315,25 +313,135 @@ int main(int argc, char** argv) { } } - //! Add InfoBlocks with statistical uncertainty from NNLOJET - if ( !STRING_ARR(InfoBlockFiles).empty() && - !STRING_ARR(InfoBlockDescr).empty() && - !STRING_ARR(InfoBlockOrders).empty()) { + //! Add InfoBlocks with statistical uncertainty from steering file + //! Default flag1=0: Statistical/numerical uncertainty + int IBFlag1 = 0; + //! Default flag2=1: Quadratic addition, alternative: 0: linear addition + int IBFlag2 = 1; + //! Default description line + std::string Default = "Please provide description!"; + if ( EXIST(InfoBlockStatUnc) ) { if ( !INT_ARR(RemoveBins).empty() ) { info["fnlo-tk-modify"]<<"Do NOT erase bins while adding InfoBlocks or vice versa! Aborted."<<endl; exit(25); } else { info["fnlo-tk-modify"]<<"Adding InfoBlocks to contributions."<<endl; } - unsigned int NFiles = STRING_ARR(InfoBlockFiles).size(); + static vector<double> dstrel_LO = DOUBLE_COL(InfoBlockStatUnc,dstrel_LO); + static vector<double> dstrel_NLO = DOUBLE_COL(InfoBlockStatUnc,dstrel_NLO); + static vector<double> dstrel_NNLO = DOUBLE_COL(InfoBlockStatUnc,dstrel_NNLO); unsigned int NDescr = STRING_ARR(InfoBlockDescr).size(); + if ( NDescr > 1 ) { + error["fnlo-tk-modify"]<<"Only one description line allowed for all blocks, aborted! NDescr = " << NDescr << endl; + exit(39); + } + int Ncontrib = table.GetNcontrib(); + int ic = 0; + for ( int i = 0; i < Ncontrib; i++ ) { + fastNLOCoeffBase* c = table.GetCoeffTable(i); + if ( fastNLOCoeffAddBase::CheckCoeffConstants(c,true) ) { + int iFlag1 = IBFlag1; + int iFlag2 = IBFlag2; + if ( EXIST(InfoBlockFlag1) ) { iFlag1 = INT(InfoBlockFlag1); } + if ( EXIST(InfoBlockFlag2) ) { iFlag2 = INT(InfoBlockFlag2); } + if ( c->IsLO() ) { + info["fnlo-tk-modify"]<<"Found LO contribution " << i << endl; + std::vector<std::string> Description; + if ( NDescr > 0 ) { + Description.push_back(STRING_ARR(InfoBlockDescr)[0]); + } else { + Description.push_back(Default); + } + if ( dstrel_LO.size() == 0 ) { + warn["fnlo-tk-modify"]<<"Found LO contribution, but no uncertainties! Nothing added." << endl; + } else if ( dstrel_LO.size() != nobs ) { + error["fnlo-tk-modify"]<<"You need the same number of uncertainties, dstrel_LO = " << dstrel_LO.size() << + ", than bins in the table, nobsbin = " << nobs << ". Aborted!" << endl; + exit(1); + } else { + c->AddCoeffInfoBlock(iFlag1,iFlag2,Description,dstrel_LO); + } + ic += 1; + } else if ( c->IsNLO() ) { + info["fnlo-tk-modify"]<<"Found NLO contribution " << i << endl; + std::vector <std:: string> Description; + if ( NDescr > 0 ) { + Description.push_back(STRING_ARR(InfoBlockDescr)[0]); + } else { + Description.push_back(Default); + } + if ( dstrel_NLO.size() == 0 ) { + warn["fnlo-tk-modify"]<<"Found NLO contribution, but no uncertainties! Nothing added." << endl; + } else if ( dstrel_NLO.size() != nobs ) { + error["fnlo-tk-modify"]<<"You need the same number of uncertainties, dstrel_NLO = " << dstrel_NLO.size() << + ", than bins in the table, nobsbin = " << nobs << ". Aborted!" << endl; + exit(1); + } else { + c->AddCoeffInfoBlock(iFlag1,iFlag2,Description,dstrel_NLO); + } + ic += 1; + } else if ( c->IsNNLO() ) { + info["fnlo-tk-modify"]<<"Found NNLO contribution " << i << endl; + std::vector <std:: string> Description; + if ( NDescr > 0 ) { + Description.push_back(STRING_ARR(InfoBlockDescr)[0]); + } else { + Description.push_back(Default); + } + if ( dstrel_NNLO.size() == 0 ) { + warn["fnlo-tk-modify"]<<"Found NNLO contribution, but no uncertainties! Nothing added." << endl; + } else if ( dstrel_NNLO.size() != nobs ) { + error["fnlo-tk-modify"]<<"You need the same number of uncertainties, dstrel_NNLO = " << dstrel_NNLO.size() << + ", than bins in the table, nobsbin = " << nobs << ". Aborted!" << endl; + exit(1); + } else { + c->AddCoeffInfoBlock(iFlag1,iFlag2,Description,dstrel_NNLO); + } + ic += 1; + } else { + info["fnlo-tk-modify"]<<"Unknown contribution " << i << endl; + info["fnlo-tk-modify"]<<"Nothing changed." << endl; + } + } + } + } + + //! Add InfoBlocks with statistical uncertainty from file (NNLOJET .dat, fnlo-tk-statunc .log, or .txt) + else if ( !STRING_ARR(InfoBlockFiles).empty() && + !STRING_ARR(InfoBlockOrders).empty() ) { + if ( !INT_ARR(RemoveBins).empty() ) { + info["fnlo-tk-modify"]<<"Do NOT erase bins while adding InfoBlocks or vice versa! Aborted."<<endl; + exit(25); + } else { + info["fnlo-tk-modify"]<<"Adding InfoBlocks to contributions."<<endl; + } + unsigned int NFiles = STRING_ARR(InfoBlockFiles).size(); + unsigned int NCols = INT_ARR(InfoBlockFileColumns).size(); unsigned int NOrders = STRING_ARR(InfoBlockOrders).size(); + unsigned int NDescr = STRING_ARR(InfoBlockDescr).size(); + if ( NFiles != NOrders ) { + error["fnlo-tk-modify"]<<"Need one order specification per file, aborted! Found NFiles = " << NFiles << ", and NOrders = " << NOrders <<endl; + exit(37); + } + unsigned int icola = 0; + unsigned int icolb = 0; + if ( NCols == 0 ) { + } else if ( NCols == 1 ) { + icola = INT_ARR(InfoBlockFileColumns)[0]; + } else if ( NCols == 2 ) { + icola = INT_ARR(InfoBlockFileColumns)[0]; + icolb = INT_ARR(InfoBlockFileColumns)[1]; + } else { + error["fnlo-tk-modify"]<<"Up to two column numbers allowed, but found more. Aborted! NCols = " << NCols <<endl; + exit(38); + } + if ( NDescr > 1 ) { + error["fnlo-tk-modify"]<<"Only one description line allowed for all blocks, aborted! NDescr = " << NDescr << endl; + exit(39); + } for ( unsigned int i = 0; i < NFiles; i++ ){ info["fnlo-tk-modify"]<<"InfoBlock file no. " << i << " is: " << STRING_ARR(InfoBlockFiles)[i] << endl; } - for ( unsigned int i = 0; i < NDescr; i++ ){ - info["fnlo-tk-modify"]<<"InfoBlock description no. " << i << " is: " << STRING_ARR(InfoBlockDescr)[i] << endl; - } for ( unsigned int i = 0; i < NOrders; i++ ){ info["fnlo-tk-modify"]<<"InfoBlock order no. " << i << " is: " << STRING_ARR(InfoBlockOrders)[i] << endl; } @@ -343,6 +451,10 @@ int main(int argc, char** argv) { for ( int i = 0; i < Ncontrib; i++ ) { fastNLOCoeffBase* c = table.GetCoeffTable(i); if ( fastNLOCoeffAddBase::CheckCoeffConstants(c,true) ) { + int iFlag1 = IBFlag1; + int iFlag2 = IBFlag2; + if ( EXIST(InfoBlockFlag1) ) { iFlag1 = INT(InfoBlockFlag1); } + if ( EXIST(InfoBlockFlag2) ) { iFlag2 = INT(InfoBlockFlag2); } if ( c->IsLO() ) { info["fnlo-tk-modify"]<<"Found LO contribution " << i << endl; int ilo = 0; @@ -354,14 +466,12 @@ int main(int argc, char** argv) { } } std::vector<std::string> Description; - if ( NDescr > 1 ) { - Description.push_back(STRING_ARR(InfoBlockDescr)[ilo]); - } else if ( NDescr > 0 ) { + if ( NDescr > 0 ) { Description.push_back(STRING_ARR(InfoBlockDescr)[0]); } else { Description.push_back(Default); } - c->AddCoeffInfoBlock(0,0,Description,STRING_ARR(InfoBlockFiles)[ilo]); + c->AddCoeffInfoBlock(iFlag1,iFlag2,Description,STRING_ARR(InfoBlockFiles)[ilo],icola,icolb); ic += 1; } else if ( c->IsNLO() ) { info["fnlo-tk-modify"]<<"Found NLO contribution " << i << endl; @@ -374,14 +484,12 @@ int main(int argc, char** argv) { } } std::vector <std:: string> Description; - if ( NDescr > 1 ) { - Description.push_back(STRING_ARR(InfoBlockDescr)[inlo]); - } else if ( NDescr > 0 ) { + if ( NDescr > 0 ) { Description.push_back(STRING_ARR(InfoBlockDescr)[0]); } else { Description.push_back(Default); } - c->AddCoeffInfoBlock(0,0,Description,STRING_ARR(InfoBlockFiles)[inlo]); + c->AddCoeffInfoBlock(iFlag1,iFlag2,Description,STRING_ARR(InfoBlockFiles)[inlo],icola,icolb); ic += 1; } else if ( c->IsNNLO() ) { info["fnlo-tk-modify"]<<"Found NNLO contribution " << i << endl; @@ -394,14 +502,12 @@ int main(int argc, char** argv) { } } std::vector <std:: string> Description; - if ( NDescr > 1 ) { - Description.push_back(STRING_ARR(InfoBlockDescr)[innlo]); - } else if ( NDescr > 0 ) { + if ( NDescr > 0 ) { Description.push_back(STRING_ARR(InfoBlockDescr)[0]); } else { Description.push_back(Default); } - c->AddCoeffInfoBlock(0,0,Description,STRING_ARR(InfoBlockFiles)[innlo]); + c->AddCoeffInfoBlock(iFlag1,iFlag2,Description,STRING_ARR(InfoBlockFiles)[innlo],icola,icolb); ic += 1; } else { info["fnlo-tk-modify"]<<"Unknown contribution " << i << endl; diff --git a/v2.3/toolkit/src/fnlo-tk-statunc.cc b/v2.3/toolkit/src/fnlo-tk-statunc.cc index 5ec74b885f35f941f901953450ba8989e2cb5e47..62d1eeccf2549291849f6c48173ed49c8bf6fd18 100644 --- a/v2.3/toolkit/src/fnlo-tk-statunc.cc +++ b/v2.3/toolkit/src/fnlo-tk-statunc.cc @@ -414,16 +414,6 @@ int main(int argc, char** argv) { } } - // //! Write out cross sections with statistical uncertainty - // yell << " " << _DSEPS << endl; - // yell << " Statistical uncertainties of the " << FixedOrderChoice << " cross section for the primary scale choice.\n"; - // yell << " " << _SSEPS << endl; - // yell << " bin cross section rel. average error of the mean\n"; - // yell << " " << _SSEPS << endl; - // for (unsigned int i=0; i<xs0.size(); i++) { - // printf("%5i %18.11e %18.11e\n",i+1,xs[i],dxs[i]/xs[i]); - // } - //! If YODA is linked, write out cross sections with statistical uncertainty in YODA format as well fastNLOLHAPDF fnlo(firsttable,PDFFile,0); @@ -439,6 +429,7 @@ int main(int argc, char** argv) { //! Get binning vector < pair < double, double > > bins = fnlo.GetObsBinsBounds(NDim-1); + //! Write out cross sections with statistical uncertainty yell << _CSEPSC << endl; shout << "fnlo-tk-statunc: Evaluating uncertainties" << endl; yell << _CSEPSC << endl;