diff --git a/LesHouches/LesHouchesFileReader.cc b/LesHouches/LesHouchesFileReader.cc --- a/LesHouches/LesHouchesFileReader.cc +++ b/LesHouches/LesHouchesFileReader.cc @@ -1,942 +1,950 @@ // -*- C++ -*- // // LesHouchesFileReader.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2019 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the LesHouchesFileReader class. // #include "LesHouchesFileReader.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include #include using namespace ThePEG; LesHouchesFileReader:: LesHouchesFileReader(const LesHouchesFileReader & x) : LesHouchesReader(x), neve(x.neve), ieve(0), LHFVersion(x.LHFVersion), outsideBlock(x.outsideBlock), headerBlock(x.headerBlock), initComments(x.initComments), initAttributes(x.initAttributes), eventComments(x.eventComments), eventAttributes(x.eventAttributes), theFileName(x.theFileName), theQNumbers(x.theQNumbers), theIncludeFxFxTags(x.theIncludeFxFxTags), theIncludeCentral(x.theIncludeCentral), theDecayer(x.theDecayer) {} LesHouchesFileReader::~LesHouchesFileReader() {} IBPtr LesHouchesFileReader::clone() const { return new_ptr(*this); } IBPtr LesHouchesFileReader::fullclone() const { return new_ptr(*this); } bool LesHouchesFileReader::preInitialize() const { return true; } void LesHouchesFileReader::doinit() { LesHouchesReader::doinit(); //cout << "theDecayer->fullName() = " << theDecayer->fullName() << endl; //if(theDecayer) { // cout << "DECAYER RETURNS TRUE" << endl; // } // are we using QNUMBERS if(!theQNumbers) return; // parse the header block and create // any new particles needed in QNUMBERS blocks string block = headerBlock; string line = ""; bool readingSLHA = false; int (*pf)(int) = tolower; unsigned int newNumber(0); do { line = StringUtils::car(block,"\r\n"); block = StringUtils::cdr(block,"\r\n"); if(line[0]=='#') continue; // are we reading the SLHA block if(readingSLHA) { // reached the end of slha block ? if(line.find(" split = StringUtils::split(line,"#"); // check for a qnumbers block transform(split[0].begin(), split[0].end(), split[0].begin(), pf); // if not contine if(split[0].find("block qnumbers")==string::npos) continue; // get name from comment string name; if(split.size()>=2) { name = StringUtils::stripws(split[1]); } else { ++newNumber; ostringstream tname; tname << "NP" << newNumber; name = tname.str(); } // extract the PDG code split = StringUtils::split(split[0]," "); istringstream is(split[2]); long PDGCode(0); is >> PDGCode; // get the charge, spin, colour and whether an antiparticle int charge(0),spin(0),colour(0),anti(0); for(unsigned int ix=0;ix<4;++ix) { line = StringUtils::car(block,"\r\n"); block = StringUtils::cdr(block,"\r\n"); int dummy[2]; istringstream is(line); is >> dummy[0] >> dummy[1]; switch (dummy[0]) { case 1: charge = dummy[1]; break; case 2: spin = dummy[1]; break; case 3: colour = dummy[1]; break; case 4: anti = dummy[1]; break; default: assert(false); } } // check if particles already exist PDPair newParticle; newParticle.first = getParticleData(PDGCode); if(newParticle.first) Throw() << "Particle with PDG code " << PDGCode << " whose creation was requested in a QNUMBERS Block" << " already exists. Retaining the original particle" << Exception::warning; if(anti) { newParticle.second = getParticleData(-PDGCode); if(newParticle.second) Throw() << "Anti-particle with PDG code " << -PDGCode << " whose creation was requested in a QNUMBERS Block" << " already exists. Retaining the original particle" << Exception::warning; if(( newParticle.first && !newParticle.second ) || ( newParticle.second && !newParticle.first ) ) Throw() << "Either particle or anti-particle with PDG code " << PDGCode << " whose creation was requested in a QNUMBERS Block" << " already exists, but not both the particle and antiparticle. " << " Something dodgy here stopping" << Exception::runerror; } // already exists continue if(newParticle.first) continue; // create the particles // particle with no anti particle if( anti == 0 ) { // construct the name if(name=="") { ostringstream temp; temp << PDGCode; name = temp.str(); } // create the ParticleData object newParticle.first = ParticleData::Create(PDGCode,name); } // particle anti-particle pair else { // construct the names string nameAnti; if(name=="") { ostringstream temp; temp << PDGCode; name = temp.str(); ostringstream temp2; temp << -PDGCode; nameAnti = temp2.str(); } else { nameAnti=name; for(string::iterator it=nameAnti.begin();it!=nameAnti.end();++it) { if(*it=='+') nameAnti.replace(it,it+1,"-"); else if(*it=='-') nameAnti.replace(it,it+1,"+"); } if(nameAnti==name) nameAnti += "bar"; } // create the ParticleData objects newParticle = ParticleData::Create(PDGCode,name,nameAnti); } // set the particle properties if(colour==1) colour = 0; newParticle.first->iColour(PDT::Colour(colour)); newParticle.first->iSpin (PDT::Spin (spin )); newParticle.first->iCharge(PDT::Charge(charge)); // register it generator()->preinitRegister(newParticle.first, "/Herwig/Particles/"+newParticle.first->PDGName()); // set the antiparticle properties if(newParticle.second) { if(colour==3||colour==6) colour *= -1; charge = -charge; newParticle.second->iColour(PDT::Colour(colour)); newParticle.second->iSpin (PDT::Spin (spin )); newParticle.second->iCharge(PDT::Charge(charge)); // register it generator()->preinitRegister(newParticle.second, "/Herwig/Particles/"+newParticle.second->PDGName()); } } // start of SLHA block ? else if(line.find("> id >> mass; // skip resetting masses on SM particles // as it can cause problems later on in event generation if(abs(id)<=6 || (abs(id)>=11 && abs(id)<=16) || abs(id)==23 || abs(id)==24) { // Throw() << "Standard model mass for PID " // << id // << " will not be changed." // << Exception::warning; block = StringUtils::cdr(block,"\r\n"); line = StringUtils::car(block,"\r\n"); continue; } // magnitude of mass for susy models mass = abs(mass); // set the mass tPDPtr particle = getParticleData(id); if(!particle) throw SetupException() << "LesHouchesFileReader::doinit() - Particle with PDG code not" << id << " not found." << Exception::runerror; const InterfaceBase * ifb = BaseRepository::FindInterface(particle, "NominalMass"); ostringstream os; os << mass; ifb->exec(*particle, "set", os.str()); // read the next line block = StringUtils::cdr(block,"\r\n"); line = StringUtils::car(block,"\r\n"); }; } // found a decay block else if(line.find("decay") == 0) { // get PGD code and width istringstream iss(line); string dummy; long parent(0); Energy width(ZERO); iss >> dummy >> parent >> iunit(width, GeV); // get the ParticleData object PDPtr inpart = getParticleData(parent); if(!inpart) { throw SetupException() << "LesHouchesFileReader::doinit() - A ParticleData object with the PDG code " << parent << " does not exist. " << Exception::runerror; return; } if ( abs(inpart->id()) == 6 || abs(inpart->id()) == 15 || abs(inpart->id()) == 23 || abs(inpart->id()) == 24 || abs(inpart->id()) == 25 ) { Throw() << "\n" "************************************************************************\n" "* Your LHE file changes the width of " << inpart->PDGName() << ".\n" "* This can cause serious problems in the event generation!\n" "************************************************************************\n" "\n" << Exception::warning; } else if (inpart->width() > ZERO && width <= ZERO) { Throw() << "\n" "************************************************************************\n" "* Your LHE file zeroes the non-zero width of " << inpart->PDGName() << ".\n" "* If " << inpart->PDGName() << " is a decaying SM particle,\n" "* this can cause serious problems in the event generation!\n" "************************************************************************\n" "\n" << Exception::warning; } // set the width inpart->width(width); if( width > ZERO ) { inpart->cTau(hbarc/width); inpart->widthCut(5.*width); inpart->stable(false); } // construct prefix for DecayModes string prefix(inpart->name() + "->"), tag(prefix),line(""); unsigned int nmode(0); // read any decay modes line = StringUtils::car(block,"\r\n"); while(line[0] != 'D' && line[0] != 'B' && line[0] != 'd' && line[0] != 'b' && line[0] != '<' && line != "") { // skip comments if(line[0] == '#') { block = StringUtils::cdr(block,"\r\n"); line = StringUtils::car(block,"\r\n"); continue; } // read decay mode and construct the tag istringstream is(line); double brat(0.); unsigned int nda(0),npr(0); is >> brat >> nda; while( true ) { long t; is >> t; if( is.fail() ) break; if( t == abs(parent) ) throw SetupException() << "An error occurred while read a decay of the " << inpart->PDGName() << ". One of its products has the same PDG code " << "as the parent particle in LesHouchesFileReader::doinit()." << " Please check the Les Houches file.\n" << Exception::runerror; tcPDPtr p = getParticleData(t); if( !p ) throw SetupException() << "LesHouchesFileReader::doinit() -" << " An unknown PDG code has been encounterd " << "while reading a decay mode. ID: " << t << Exception::runerror; ++npr; tag += p->name() + ","; } if( npr != nda ) throw SetupException() << "LesHouchesFileReader::doinit() - While reading a decay of the " << inpart->PDGName() << " from an SLHA file, an inconsistency " << "between the number of decay products and the value in " << "the 'NDA' column was found. Please check if the spectrum " << "file is correct.\n" << Exception::warning; // create the DecayMode if( npr > 1 ) { if( nmode==0 ) { generator()->preinitInterface(inpart, "VariableRatio" , "set","false"); if(inpart->massGenerator()) { ok = false; Throw() << inpart->PDGName() << " already has a MassGenerator set" << " this is incompatible with using QNUMBERS " << "Use\n" << "set " << inpart->fullName() << ":Mass_generator NULL\n" << "to fix this." << Exception::warning; } if(inpart->widthGenerator()) { ok = false; Throw() << inpart->PDGName() << " already has a WidthGenerator set" << " this is incompatible with using QNUMBERS " << "Use\n" << "set " << inpart->fullName() << ":Width_generator NULL\n" << "to fix this." << Exception::warning; } unsigned int ntemp=0; for(DecaySet::const_iterator dit = inpart->decayModes().begin(); dit != inpart->decayModes().end(); ++dit ) { if((**dit).on()) ++ntemp; } if(ntemp!=0) { ok = false; Throw() << inpart->PDGName() << " already has DecayModes" << " this is incompatible with using QNUMBERS " << "Use\n" << "do " << inpart->fullName() << ":SelectDecayModes none\n" << " to fix this." << Exception::warning; } } inpart->stable(false); tag.replace(tag.size() - 1, 1, ";"); DMPtr dm = generator()->findDecayMode(tag); if(!theDecayer) Throw() << "LesHouchesFileReader::doinit() Decayer must be set using the " << "LesHouchesFileReader:Decayer" << " must be set to allow the creation of new" << " decay modes." << Exception::runerror; if(!dm) { dm = generator()->preinitCreateDecayMode(tag); if(!dm) Throw() << "LesHouchesFileReader::doinit() - Needed to create " << "new decaymode but one could not be created for the tag " << tag << Exception::warning; } generator()->preinitInterface(dm, "Decayer", "set", theDecayer->fullName()); ostringstream br; br << setprecision(13) << brat; generator()->preinitInterface(dm, "BranchingRatio", "set", br.str()); generator()->preinitInterface(dm, "Active", "set", "Yes"); if(dm->CC()) { generator()->preinitInterface(dm->CC(), "BranchingRatio", "set", br.str()); generator()->preinitInterface(dm->CC(), "Active", "set", "Yes"); } ++nmode; } tag=prefix; // read the next line block = StringUtils::cdr(block,"\r\n"); line = StringUtils::car(block,"\r\n"); }; if(nmode>0) { inpart->update(); if(inpart->CC()) inpart->CC()->update(); } } } // start of SLHA block ? else if(line.find("() << "The file associated with '" << name() << "' does not contain a " << "proper formatted Les Houches event file. The events may not be " << "properly sampled." << Exception::warning; } //vector LesHouchesFileReader::optWeightNamesFunc() { return optionalWeightsNames; } vector LesHouchesFileReader::optWeightsNamesFunc() { return optionalWeightsNames; } void LesHouchesFileReader::open() { if ( filename().empty() ) throw LesHouchesFileError() << "No Les Houches file name. " << "Use 'set " << name() << ":FileName'." << Exception::runerror; cfile.open(filename()); if ( !cfile ) throw LesHouchesFileError() << "The LesHouchesFileReader '" << name() << "' could not open the " << "event file called '" << theFileName << "'." << Exception::runerror; cfile.readline(); if ( !cfile.find(" attributes = StringUtils::xmlAttributes("LesHouchesEvents", cfile.getline()); LHFVersion = attributes["version"]; if ( LHFVersion.empty() ) return; bool readingHeader = false; bool readingInit = false; headerBlock = ""; // char (cwgtinfo_weights_info[250][15]); string hs; // int cwgtinfo_nn(0); // Loop over all lines until we hit the tag. bool readingInitWeights = false, readingInitWeights_sc = false; string weightinfo; while ( cfile.readline() ) { // found the init block for multiple weights if(cfile.find(""; erase_substr(sub, str_arrow); scalename = sub; break; } ++ws; } while (isc); /* now get the relevant information * e.g. scales or PDF sets used */ string startDEL = "'>"; //starting delimiter string stopDEL = ""; //end delimiter int firstLim = hs.find(startDEL); //find start of delimiter if(firstLim == -1) { startDEL = ">"; firstLim = hs.find(startDEL); } // unsigned lastLim = hs.find(stopDEL); //find end of delimitr string scinfo = hs.substr(firstLim); //define the information for the scale erase_substr(scinfo,stopDEL); erase_substr(scinfo,startDEL); scinfo = StringUtils::stripws(scinfo); /* fill in the map * indicating the information to be appended to each scale * i.e. scinfo for each scalname */ scalemap[scalename] = scinfo.c_str(); string str_id = "id="; string str_prime = "'"; erase_substr(scalename, str_id); erase_substr(scalename, str_prime); optionalWeightsNames.push_back(scalename); } } if ( cfile.find("") ) { //cout << "found init block" << endl; // We have hit the init block, so we should expect to find the // standard information in the following. But first check for // attributes. initAttributes = StringUtils::xmlAttributes("init", cfile.getline()); readingInit = true; cfile.readline(); if ( !( cfile >> heprup.IDBMUP.first >> heprup.IDBMUP.second >> heprup.EBMUP.first >> heprup.EBMUP.second >> heprup.PDFGUP.first >> heprup.PDFGUP.second >> heprup.PDFSUP.first >> heprup.PDFSUP.second >> heprup.IDWTUP >> heprup.NPRUP ) ) { heprup.NPRUP = -42; LHFVersion = ""; return; } heprup.resize(); for ( int i = 0; i < heprup.NPRUP; ++i ) { cfile.readline(); if ( !( cfile >> heprup.XSECUP[i] >> heprup.XERRUP[i] >> heprup.XMAXUP[i] >> heprup.LPRUP[i] ) ) { heprup.NPRUP = -42; LHFVersion = ""; return; } } } if ( cfile.find("> sub; if(we==2) { npLO = atoi(sub.c_str()); } if(we==5) { npNLO = atoi(sub.c_str()); } ++we; } while (ievat); optionalnpLO = npLO; optionalnpNLO = npNLO; std::stringstream npstringstream; npstringstream << "np " << npLO << " " << npNLO; std::string npstrings = npstringstream.str(); /* the FxFx merging information * becomes part of the optionalWeights, labelled -999 * for future reference */ if(theIncludeFxFxTags) optionalWeights[npstrings.c_str()] = -999; if ( !cfile.readline() ) return false; // The first line determines how many subsequent particle lines we // have. if ( !( cfile >> hepeup.NUP >> hepeup.IDPRUP >> hepeup.XWGTUP >> hepeup.SCALUP >> hepeup.AQEDUP >> hepeup.AQCDUP ) ) return false; hepeup.resize(); // Read all particle lines. for ( int i = 0; i < hepeup.NUP; ++i ) { if ( !cfile.readline() ) return false; if ( !( cfile >> hepeup.IDUP[i] >> hepeup.ISTUP[i] >> hepeup.MOTHUP[i].first >> hepeup.MOTHUP[i].second >> hepeup.ICOLUP[i].first >> hepeup.ICOLUP[i].second >> hepeup.PUP[i][0] >> hepeup.PUP[i][1] >> hepeup.PUP[i][2] >> hepeup.PUP[i][3] >> hepeup.PUP[i][4] >> hepeup.VTIMUP[i] >> hepeup.SPINUP[i] ) ) return false; //print momenta for debugging // cout << hepeup.PUP[i][0] << " " << hepeup.PUP[i][1] << " " << hepeup.PUP[i][2] << " " << hepeup.PUP[i][3] << " " << hepeup.PUP[i][4] << endl; if(std::isnan(hepeup.PUP[i][0])||std::isnan(hepeup.PUP[i][1])|| std::isnan(hepeup.PUP[i][2])||std::isnan(hepeup.PUP[i][3])|| std::isnan(hepeup.PUP[i][4])) throw Exception() << "nan's as momenta in Les Houches file " << Exception::eventerror; if(hepeup.MOTHUP[i].first -1==i || hepeup.MOTHUP[i].second-1==i) { throw Exception() << "Particle has itself as a mother in Les Houches " << "file, this is not allowed\n" << Exception::eventerror; } } // Now read any additional comments and named weights. // read until the end of rwgt is found bool readingWeights = false, readingaMCFast = false, readingMG5ClusInfo = false; while ( cfile.readline() && !cfile.find("")) { if(cfile.find("> sub; if(wi==1) { string str_arrow = ">" ; erase_substr(sub, str_arrow); weightName = sub; } if(wi==2) weightValue = atof(sub.c_str()); ++wi; } while (iss); // store the optional weights found in the temporary map optionalWeightsTemp[weightName] = weightValue; } /* reading of aMCFast weights */ if(readingaMCFast) { std::stringstream amcfstringstream; amcfstringstream << "aMCFast " << cfile.getline(); std::string amcfstrings = amcfstringstream.str(); string str_newline = "\n"; erase_substr(amcfstrings,str_newline); optionalWeights[amcfstrings.c_str()] = -111; //for the aMCFast we give them a weight -111 for future identification } /* read additional MG5 Clustering information * used in LO merging */ if(readingMG5ClusInfo) { string hs = cfile.getline(); string startDEL = ""; //end delimiter unsigned firstLim = hs.find(startDEL); //find start of delimiter // unsigned lastLim = hs.find(stopDEL); //find end of delimitr string mg5clusinfo = hs.substr(firstLim); //define the information for the scale erase_substr(mg5clusinfo,stopDEL); erase_substr(mg5clusinfo,startDEL); string str_arrow = ">"; erase_substr(mg5clusinfo,str_arrow); string str_quotation = "\""; erase_substr(mg5clusinfo,str_quotation); optionalWeights[mg5clusinfo.c_str()] = -222; //for the mg5 scale info weights we give them a weight -222 for future identification } //store MG5 clustering information if(cfile.find("::const_iterator it=optionalWeightsTemp.begin(); it!=optionalWeightsTemp.end(); ++it){ //to avoid issues with inconsistencies of defining the weight ids, remove "" and '' string id_1 = it->first; erase_substr(id_1, str_quote); erase_substr(id_1, str_doublequote); for (map::const_iterator it2=scalemap.begin(); it2!=scalemap.end(); ++it2){ //find the scale id in the scale information and add this information string id_2 = it2->first; erase_substr(id_2, str_quote); erase_substr(id_2, str_doublequote); if(id_1==id_2) { string info = it2->second; string str_newline = "\n"; erase_substr(info, str_newline); //set the optional weights optionalWeights[info] = it->second; } } } /* additionally, we set the "central" scale * this is actually the default event weight */ string central = "central"; if (theIncludeCentral) optionalWeights[central] = hepeup.XWGTUP; if ( !cfile ) return false; return true; } void LesHouchesFileReader::close() { cfile.close(); } void LesHouchesFileReader::persistentOutput(PersistentOStream & os) const { os << neve << LHFVersion << outsideBlock << headerBlock << initComments << initAttributes << eventComments << eventAttributes << theFileName << theQNumbers << theIncludeFxFxTags << theIncludeCentral << theDecayer ; } void LesHouchesFileReader::persistentInput(PersistentIStream & is, int) { is >> neve >> LHFVersion >> outsideBlock >> headerBlock >> initComments >> initAttributes >> eventComments >> eventAttributes >> theFileName >> theQNumbers >> theIncludeFxFxTags >> theIncludeCentral >> theDecayer; ieve = 0; } ClassDescription LesHouchesFileReader::initLesHouchesFileReader; // Definition of the static class description member. void LesHouchesFileReader::Init() { static ClassDocumentation documentation ("ThePEG::LesHouchesFileReader is an base class to be used for objects " "which reads event files from matrix element generators. This class is " "able to read plain event files conforming to the Les Houches Event File " "accord."); static Parameter interfaceFileName ("FileName", "The name of a file containing events conforming to the Les Houches " "protocol to be read into ThePEG. A file name ending in " ".gz will be read from a pipe which uses " "zcat. If a file name ends in | the " "preceeding string is interpreted as a command, the output of which " "will be read through a pipe.", &LesHouchesFileReader::theFileName, "", false, false); interfaceFileName.fileType(); interfaceFileName.rank(11); static Switch interfaceQNumbers ("QNumbers", "Whether or not to read search for and read a QNUMBERS" " block in the header of the file.", &LesHouchesFileReader::theQNumbers, false, false, false); static SwitchOption interfaceQNumbersYes (interfaceQNumbers, "Yes", "Use QNUMBERS", true); static SwitchOption interfaceQNumbersNo (interfaceQNumbers, "No", "Don't use QNUMBERS", false); static Switch interfaceIncludeFxFxTags ("IncludeFxFxTags", "Include FxFx tags", &LesHouchesFileReader::theIncludeFxFxTags, false, true, false); static SwitchOption interfaceIncludeFxFxTagsYes (interfaceIncludeFxFxTags, "Yes", "Use the FxFx tags", true); static SwitchOption interfaceIncludeFxFxTagsNo (interfaceIncludeFxFxTags, "No", "Don't use the FxFx tags", false); static Switch interfaceIncludeCentral ("IncludeCentral", "Include definition of central weight", &LesHouchesFileReader::theIncludeCentral, false, true, false); static SwitchOption interfaceIncludeCentralYes (interfaceIncludeCentral, "Yes", "include definition of central weight", true); static SwitchOption interfaceIncludeCentralNo (interfaceIncludeCentral, "No", "Don't include definition of central weight", false); static Reference interfaceDecayer ("Decayer", "Decayer to use for any decays read from the QNUMBERS Blocks", &LesHouchesFileReader::theDecayer, false, false, true, true, false); } void LesHouchesFileReader::erase_substr(std::string& subject, const std::string& search) { size_t pos = 0; while((pos = subject.find(search, pos)) != std::string::npos) { subject.erase( pos, search.length() ); } } diff --git a/NEWS b/NEWS --- a/NEWS +++ b/NEWS @@ -1,1009 +1,1015 @@ ThePEG News -*- outline -*- =============================== Numbered bugs can be found at https://phab.hepforge.org/TNN where NN is the bug number. The latest version of ThePEG can be found at http://www.thep.lu.se/ThePEG or at https:///herwig.hepforge.org/ +* ThePEG-2.2.2 release: 2021-01-22 + +** add meta mechanism to event record + +** add option of preloading yoda files and remove support for rivet-1 and aida + * ThePEG-2.2.1 release: 2019-04-10 ** new release solely to keep Herwig and ThePEG version numbers in sync * ThePEG-2.2.0 release: 2019-12-12 ** Improvements to the unit templates, building with gcc 4.8 is no longer possible ** Several smaller bug fixes and additions, to allow new features in Herwig 7.2 * ThePEG-2.1.7 release: 2020-04-10 ** new release solely to keep Herwig and ThePEG version numbers in sync * ThePEG-2.1.6 release: 2019-12-11 ** update boost.m4 to allow compliation with new boost versions ** increase the eps parameter for boosts of wavefunctions ** fixes for gcc4.8 and C++17 thanks to Ivan Razumov (GENSER) ** Added inelasticity cut in SimpleDISCut (thanks Andrii Verbytskyi ) * ThePEG-2.1.5 release: 2019-04-04 ** Improvements to template instantation of templates with gcc9 and icc T23 ** Change in assignment operator definition, add delete, to avoid warnings with gcc9 ** Minor changes for compilation with gcc8,9, icc and clang * ThePEG-2.1.4 release: 2018-28-06 ** Improvements to helicity libraries to support more BSM models ** Added FixedTargetLuminosity for fixed target collisions * ThePEG-2.1.3 release: 2018-04-05 ** Use std::array<> where suitable * ThePEG-2.1.2 release: 2017-11-01 ** Allow LHAPDF6 interface to return photon pdf. ** Sign fix in GeneralVVSVertex ** Fix multiple weight reading from LHE files ** Preparations for Rivet 3 * ThePEG-2.1.1 release: 2017-07-14 ** Added missing evaluate() option to GeneralVVSVertex ** More robust reading of LHE files ** Warn about duplicate PDG names in user input ** Write current timestamp into log files ** Calling 'read' in an input file will no longer change the repository directory you're in * ThePEG-2.1.0 release: 2017-05-19 ** Transition to C++-11 code, building with C++-98 is no longer possible ** Interfaces with two choices now consistently accept Yes/No as answer ** Several smaller bug fixes and additions, to allow new features in Herwig 7.1 * ThePEG-2.0.4 release: 2016-10-25 ** Default weight explicitly labelled The nominal event weight is now labelled "Default" and is always the first weight to be inserted into HepMC. (Note that HepMC 2.06.09 may not preserve this ordering.) ** LesHouches event weights handling There are now two options for the treatment of optional weights coming in from a Les Houches file, available through "LesHouchesEventHandler:WeightNormalization". Either they are normalized to the "Default" weight ("Normalized"), or they are given as cross-sections in pb ("CrossSection"). The default behaviour is "Normalized". ** JetFinder in Cuts objects writes information The selected JetFinder is explicitly listed in the debug output * ThePEG-2.0.3 release: 2016-07-21 ** LesHouches reader Removed FxFx-specific weights from default output ** FuzzyTheta cut Fixed missing division by width ** NaN check for incoming momenta from ME providers ** Build system Compatibility with gcc-6, improved configure macros. * ThePEG-2.0.2 release: 2016-04-27 ** LesHouches reader Fixed cross-section issue for LH event files using weight scheme 3 * ThePEG-2.0.1 release: 2016-02-17 ** Electroweak scheme fixes Schemes 2 and 5 now set GF, scheme 7 calculates mw correctly. ** New ParVector 'clear' command ** More visible version information The version number of ThePEG is written into the log file. ** Rivet exception handling Any exceptions from Rivet/YODA are converted into warning exceptions. This can invalidate plots, please check the logfile. ** Rapidity limits fixed in FuzzyTheta cuts ** Doxygen updated to version 1.8 ** Java interface fixed * ThePEG-2.0.0 release: 2015-12-04 ** Improvements for NLO A number of improvements and bug fixes have been made in the course of (Herwig 7) NLO development. ** Improvements for spin correlations A number of improvements have been made to handle spin correlations in the parton shower. ** Setup file mechanism An isolated event generator can now be modified prior to generating events by providing an additional input file; also support for re-initalizing isolated event generators has been added in this course to support Herwig 7's parallel integration modes. ** Handling of weighted LHE events The handling of LHE events with variable weight (codes 3 and 4) has now been enabled, along with the full handling of multiple event weights down to the HepMC output. ** LHAPDF6 handling A problem in initializing LHAPDF6 has been fixed. ** Heavy ion collisions A dedicated HepMC interface for heavy ion collisions has been introduced. ** First steps for unit testing Unit testing has been introduced for a number of components based on boost's unittesting framework. They are only enabled, when configured --with-boost . * ThePEG-1.9.2 release: 2014-07-07 ** Better support for LHAPDF 6 The LHAPDF interface now determines if version 6 is available and makes use of new LHAPDFv6 features. * ThePEG-1.9.1 release: 2014-04-30 ** Build fix for SLC6 Fixed problems with a missing header in ACDCGen.h that affected some builds on SLC6 ** Build fix for Rivet 2.1.1 Adapted to Rivet 2.1.1 changed header file layout ** OS X Mavericks build fix for readline The problem of infinite hangs when linking against libreadline in Mavericks has been resolved. * ThePEG-1.9.0 release: 2013-10-28 ** Rapidity edge cases Changed behavior for vanishing pt in eta() and vanishing mt in rapidity(): no error is thrown, but a ridculously large rapidity is returned instead. ** Command line interface Adding a tag of the form "#first-last" when running a MultiEventGenerator now allows to run a subset of the runs ** Pre- and Posthandlers A warning is now issued if the same handler is inserted twice into a pre- or posthandler list. ** Statistics fixes Various issues for the handling of statistics in the presence of negative weights have been fixed. ** Rivet search paths A new interface RivetAnalysis:Paths allows more search paths to be added to the analysis finding routine. This supplements the RIVET_ANALYSIS_PATH environment variable. ** Vertex classes A chain of sign and prefactor fixes has been implemented in the Vertex classes to correct various vertex structures in Herwig++ conversions of Feynrules models. ** Custom XComb objects Matrix element classes can built with customized XComb objects. ** Reweighting of SubProcessGroups and projections Matrix element groups can reweight contributions in a SubProcessGroup and select one of the dependent subprocesses to be used as the hard one instead of the group in total. ** Fuzzy cuts Jet cuts in the jet cut framework can be assigned a fuzzy shape to improve the stability of NLO calculations. ** NLORivetAnalysis, NLOHepMCFile The NLOHepMCFile and NLORivetAnalysis classes have been added to allow Rivet to analyse NLO calculations which perform output for the subtracted real emission matrix elements as ThePEG::SubProcessGroups; individual subprocesses in these groups are flagged as correlated by being assigned the same event number. ** C++-11 testing To help with the coming transition to C++-11, we provide the new --enable-stdcxx11 configure flag. Please try to test builds with this flag enabled and let us know any problems, but do not use this in production code yet. In future releases, this flag will be on by default. * ThePEG-1.8.3 release: 2013-02-22 ** Fixed HepMC status code assignment The incoming line to a technical vertex that only changes a particle's momentum is now labelled 11 (MC-internal) instead of 2 (unstable). * ThePEG-1.8.2 release: 2013-01-30 ** Changes to scheduled event dumps Set the EventGenerator option "KeepAllDumps" to keep all dump files of a run, labelled by event number. These scheduled dumps now produce a cleaned generator state between events. ** SLHA block bug fix In some circumstances, the end of an SLHA block was not parsed correctly. * ThePEG-1.8.1 release: 2012-10-15 ** Repository changes *** Stable flag The behaviour of the 'Particle:Stable' flag has changed slightly. When all decaymodes are removed from a particle, the flag will be automatically set to stable. However, it still needs to be set unstable by hand when new decaymodes are added to a formerly stable particle. *** Search paths The search paths for reading input files are now saved in the repository. ** Rivet analysis YODA capable The configure time check for Rivet will now detect the version to check if it is YODA capable. ** NLOHepMCFile An AnalysisHandler, similar to the standard HepMCFile has been added to allow analysing plain NLO calculations by communicating each member of a subprocess group (typically real emission and subtraction contributions) to a HepMC file; until Rivet supports the correct treatment of errors in this case, the subprocess group members are written out as events independent of each other. ** Jet cuts A more flexible framework for cuts on jets has been added, including jet finders (either builtin or via an optional link to fastjet). All kinds of exclusive and inclusive jet cuts, as well as jet vetoes on the level of the hard process are now supported. ** Tree2toNDiagram A bug has been fixed in Tree2toNDiagram::isSame(tcDiagPtr, map&) method to take into account symmetries when swapping two external legs attached to the same vertex ** Quark thresholds in alphaS Support for setting quark masses for threshold matching in couplings deriving from AlphaSBase independently of masses used in ParticleData objects has been added. The O1AlphaS implementation is fully aware of this enhancement. ** LesHouches reader fixes The LesHouches readers have been improved, and will spot many more consistency errors, such as coloured Standard Model leptons, or helicity values outside of [-1,1] and 9. ** Spinor enhancements The LorentzSpinor class has gained member functions for projection operations and the contraction with the sigma-mu-nu commutator. * ThePEG-1.8.0 release: 2012-05-21 ** NLO support *** ThePEG now includes structures to ease implementing next-to-leading order (NLO) calculations as well as for interfacing external matrix element codes through runtime interfaces. Particularly, the newly introduced MEGroup and accompanying StdXCombGroup, StdDependentXComb and SubProcessGroup classes provide the functionality required by subtraction approaches to higher orders. A general interface for cutting on reconstructed jets as required by higher-order calculations is included, along with an implementation of kt, Cambridge-Aachen and anti-kt jet finding as relevant for NLO calculations. Hard process implementations deriving from MEBase are no longer limited to the evaluation of PDFs by PartonExtractor objects, thus allowing for a more flexible and more stable implementation of finite collinear contributioins appearing in the context of higher order corrections. *** The generation of phasespace points for the hard subprocess has been made more flexible, particularly to allow generation of incoming parton momenta by the hard matrix element as is typically done by phasespace generators provided with fixed-order codes. Along with this change, generation of the phasespace point does not need to take place in the centre-of-mass system of the incoming partons. *** Various helpers have been added to MEBase and dependent code along with the improvements described above, including simple functionality required for caching intermediate results. *** Tree2toNDiagram supports merging of two external legs, yielding another Tree2toNDiagram object to assist in determining subtraction terms required for a particular process ** Support for SU(3)-sextet colour lines ** Named weights support, also in HepMC Named, optional weights on top of the usual event weight are now fully supported; this includes their communication to HepMC events as well as their parsing from the extented Les Houches file format drafted at the Les Houches workshop 2009. ** Complex masses supported in Helicity code ** Several minor fixes and enhancements * ThePEG-1.7.3 release: 2012-03-05 The only changes are in LesHouches.so, now at version 14. ** Spin information Spin correlation information will now be set up correctly for tau leptons. To go back to the old behaviour, set LesHouchesReader:IncludeSpin No ** Consistency checks Catch broken input where mother-daughter relations are circular. * ThePEG-1.7.2 release: 2011-11-01 ** HepMC configuration Clarified at configure time that HepMC versions before 2.05 are not officially supported. ** Rivet configuration Rivet builds with external header dependencies are now correctly recognized. ** Helicity vertex consistency To help debugging, the addToList() function for registering particle lines connected to a vertex now checks for electric charge conservation at the vertex. Additionally, the specified QED/QCD order of the interaction is checked. ** Ticket #355: Global library list for resume functionality The list of global libraries is now correctly included in the dump file, to fix the functionality of resuming a run from regular checkpoints. * ThePEG-1.7.1 release: 2011-06-20 ** Ticket #238: Self-consistent electroweak schemes The default behaviour of Herwig++ is to use the best values for all the electroweak parameters, which can be inconsistent. Optionally now, a self-consistent electroweak scheme can be chosen in the parameter Model:EW/Scheme. Note that values of 'EW/Scheme' away from the default have not received the same amount of testing. ** Ticket #338: Fixed reporting of floating point exceptions This was an issue whenever LHAPDF libraries were linked in. ** Fixed cTau() behaviour The new behavior is that if both width and lifetime are zero, cTau() returns zero for unstable particles and MaxLength for stable ones. ** MaxErrors The cutoff can be disabled by setting MaxErrors to -1. ** StdOut redirect CurrentGenerator::Redirect now does not redirect to to the internal stream in EventGenerator if the useStdout flag has been set. ** Run name tags Aded possibility to add a tag to the run name when running with the '-t' option. One run file can thus be run with different seeds and result in different output files. ** Exception names EventGenerator tries to convert exception type names into human-readable form. Currently this only works for gcc-compatible compilers. ** Repository API changes Instead of printing an error message to cerr, the Repository::load() and Repository::read(filename,os) commands now behave like the other repo commands and return an error string. This allows --exitonerror to work correctly for load() and read(). Users of these functions need to send the string to cerr themselves if the old output behaviour is required. Repository::read(is, os, prompt) is unchanged. ** HepMC precision The precision() option for HepMC GenEvent is now available as an interface in the HepMCFile analysis handler. ** gcc-4.6 The build has now also been tested with gcc 4.6.0. * ThePEG-1.7.0 release: 2011-02-08 ** Behaviour *** Cross-section information The .out file contains a better estimate of the cross-section directly from the phase space sampler. It should be reliable if not too many events were vetoed during the later phases of the production. *** Simpler decay mode selection Instead of having to turn off every mode individually, the new 'SelectDecayModes' interface allows commands such as tbar:SelectDecayModes none tbar:SelectDecayModes tbar->nu_ebar,e-,bbar; tbar->nu_mubar,mu-,bbar; Use 'PrintDecayModes' to list the available choices. *** Rivet interface The interface will check before the run if all chosen analyses are actually available in Rivet, and will only call finalize() if any events have been generated. Event weights are now passed correctly into Rivet. *** Les Houches QNUMBER support QNUMBER particle creation support has been added to the Les Houches reader. *** Debug level If a debug level is set on the command line, it will always be used. *** Abort preserves events A hard abort exception during event generation will try to finalize as best as it can, thus preserving information about the run until this point. *** Progress log There is a new ThePEG::ProgressLog analysis handler, which prints timing information about the run to screen. *** Graphviz The event graph can show missing momentum information. *** CRLF line endings (Windows-style) File readers can now cope with files that have CRLF (Windows) and CR (Mac) line endings. ** Structure *** Diffraction support To provide support for the simulation of diffractive events, the LeptonLeptonRemnant class was renamed to UnresolvedRemnant, the WeizsackerWilliams PDF extended, BudnevPDF added, and PartonExtractor modified to allow separate override PDFs for each beam. *** Polarized beams The PolarizedBeamParticleData class permits a spin polarization choice on the incoming particles. *** Mixing particles The MixedParticleData class supports mixed particle states. *** SpinBase The SpinBase class has been removed, SpinInfo is now the base class. *** Vertex classes Several new vertices were added for BSM physics, and some minor bugs fixed with existing classes. *** PID type Particle IDs now have their own type, 'PID'. It can only be converted to 'long' to avoid unsigned int overflow errors on 64bit machines. *** ClassDescription There is a simpler way to register ClassDescription information with the Repository: the DescribeClass<> template. No information is needed in the class headers anymore, reducing build dependencies. New code will use this method now, older classes will be migrated in future. *** LWH histogramming The built-in implementation of AIDA histogramming has been restructured slightly to decouple the implementation from the interface. ** Build system *** Silent build rules 'make' now builds silently, to improve readability. To get the old behaviour, run 'make V=1' *** zlib Integrated zlib to read compressed files, such as Les Houches events. *** Libtool 2.4 Will fix some issues on OS X. *** gsl check This check now also works with Rivet ** Bug fixes *** Ticket #286: Implementation of Cuts in LesHouchesReader Les Houches events are given in the lab frame, but QCDCuts expects to be given momenta in the parton-parton cmf. Boost added. *** Ticket #303: AlphaS thresholds AlphaS thresholds were fixed for exactly massless quarks. *** Ticket #304: Length units in HepMC *** Ticket #309: Particle initialization order *** Ticket #310: LeptonLeptonPDF fix * ThePEG-1.6.1 release: 2009-12-11 ** Ticket #292: Spin correlation fix to stabilize tau decay numerics We have restructured the spin correlation code to stabilize the numerical problems seen in tau decays. ** Speed increase Improved decay chain handling speeds up a typical run by a few percent. ** Exception logging The log file lists different exception types individually again instead of grouping them all by severity only. This was broken in the last few releases. * ThePEG-1.6.0 release: 2009-11-19 ** Helicity amplitudes The main change in this release is a streamlining of the helicity amplitude code. If you have self-written vertex classes, you will need to adapt them slightly: *** Spinor representation All spinors are now in HELAS representation, the optional representation switching has been removed. *** Adding particles to the vertex Instead of calling setList(v1, v2, v3) with three vectors of particles that need to be filled in sync, call addToList(p1, p2, p3) repeatedly. *** Name changes To make them consistent with the rest of ThePEG, the names of functions starting with get...() have changed: getFoo() becomes foo(); setFoo(...) becomes foo(...). ** Vector code changes *** Vector3 has been renamed ThreeVector This makes it consistent with the other vector classes. *** LorentzVector::mag()/mag2() removed, to reduce confusion Use the equivalent ::m()/m2() instead to get E^2-p^2. LorentzVector::rho2() == ThreeVector::mag2() == p^2 *** Calculating with a zero-length vector Previously, mathematically undefined values were arbitrarily set to 0 in this situation. Now an assertion failure is triggered in debug mode. Only the azimuthal angle phi() still returns 0. ** Environment variable interpretation removed The programs setupThePEG and runThePEG are not wrapped in shell scripts anymore. All usage of environment variables at runtime has been removed. To influence the behaviour of ThePEG, you will need to use explicit command line flags, or calls to Repository:: functions. ** StandardModelBase now provides G_Fermi InvEnergy2 StandardModelBase::fermiConstant() const; provides the PDG 2006 value of G_Fermi. ** Allow particle width cuts to be unset Setting a negative value for the lower or upper width cut removes that bound on the width. ** Ticket #271: Decay mode names are normalized Previously, the ordering of decay products in a decay mode specifier needed to match ThePEG's internal ordering exactly. This is now done automatically. ** Ticket #273: NoPDF and LesHouches reader The LesHouches reader has been fixed for the case where no PDFs are used in the LHE file. ** Ticket #275: doinitrun() ordering The ACDC sampler now ensures that initrun() of all ME objects is run first. ** Ticket #277: Repository command help texts ThePEG's repository command language now includes a 'help' functionality. ** Redirection of various files The log file stream now goes to stdout instead of stderr when EventGenerator:UseStdout is set. ** Readline support can be disabled Using the configure switch '--disable-readline', the linking of libreadline can be suppressed. This can be useful for batch code that will never be used interactively, if the number of linked-in libraries is a problem. ** --with-LHAPDF configure flag Previously, the full LHAPDF path including lib/ needed to be specified. This now also works with the more common usage of just LHAPDF's prefix path. ** Rivet analysis plugin The Rivet analysis output files are now named consistently like the other output from a run, with the run name as prefix. ** Graphviz event visualization This is now independent of HepMC. The call to void ThePEG::printGraphviz(ostream & os, tcEventPtr event); on any event will output a Graphviz file to the stream, suitable for interpretation with the 'dot' tool. It shows a visualization of the generator's internal event structure, useful for debugging. This is an initial version, feedback is welcome! ** Fixed compatibility with older HepMC versions A problem with rejecting a missing HepMC unit implementation was fixed. * ThePEG-1.5.0 release: 2009-09-01 ** New ZERO object The ZERO object can be used to set any dimensionful quantity to zero. This avoids explicit constructs like 0.0*GeV. ** Readline interface The interactive repository access now uses the readline library, providing a command history and easier command line editing. ** Syntax highlighting We now have a syntax highlighting mode for emacs. To enable it, use 'M-x ThePEG-repository-mode' on any .in file. ** Configure information Important configuration information is listed at the end of the 'configure' run and in the file 'config.thepeg'. Please provide this file in any bug reports. ** Rivet interface ThePEG now supports Rivet internally as an AnalysisHandler if it's available. ** HepMC tools; CLHEP support removed The HepMC file output and graphviz event view have migrated from Herwig++ to ThePEG. The deprecated CLHEP support has been removed. ** Exception specifiers removed Client code changes are needed in doinit() etc. Simply remove the exception specifier after the function name. ** Support for HepMC 2.05 *** New features ThePEG now supports cross-section output, PDF information and unit specifications if they are available in HepMC. *** IO_ASCII Support for the deprecated IO_ASCII format has been removed. *** Status codes ThePEG uses codes 1, 2 and 11 according to the HepMC agreement. ** Redirection of .out, .log and .tex to stdout Set 'EventGenerator:UseStdout' to 'Yes' and (almost) all output will be streamed to stdout instead of files. ** Ticket #248: Les Houches reader The cross-section information is now reported correctly when reading several files. ** Cuts output If the debug level is set > 0, the current set of cuts is prepended to the logfile. ** Preweighting of matrix elements A segfault when using preweights was fixed. Preweights are now correctly included in handler statistics. ** Other technical changes *** Colour line improvements *** PDFsets.index search improved *** Ticket #232: Java check on OS X now works headless *** Running couplings restructured *** LeptonLeptonRemnant iprovements to support GammaGamma *** WaveFunction constructors streamlined *** VertexBase now provides sin2ThetaW etc. * ThePEG-1.4.2 release: 2009-05-08 ** Ticket #242 Fixed a compiler problem that showed up on openSUSE 10.2, g++ 4.1.2. A source line was omitted if the optimization level was higher than -O1. ** User interaction Dump file generation can now be disabled completely by setting EventGenerator:DumpPeriod to -1. * ThePEG-1.4.1 release: 2009-03-31 ** User interaction Error messages have been clarified in BaseRepository and StandardEventHandler ** Les Houches files File readers are more robust, with clearer messages when things go wrong. ** Floating point issues fixed in ThreeVector and VertexBase. ** HepMC converter Fixed PDF choice for asymmetric beams. ** Libtool Updated to version 2.2.6 * ThePEG-1.4.0 release: 2008-12-02 ** Efficiency improvements The LorentzVector class, the helicity amplitude code and PDF lookups have been profiled and restructured to eliminate speed bottlenecks. ** Deep inelastic scattering Support for DIS is now implemented in ThePEG. ** New rapidity cut Added alternative pT cut that cuts on rapidity rather than pseudorapidty as the existing one fails for zero-pt massive particles. ** HepMC units support Users of HepMC versions > 2.04 get full units support. The HepMC GenEvent object has the correct units set. ** Support for HepMC PdfInfo data Users of the HepMC converter can now fill the PdfInfo data, which has been available since HepMC 1.28.00. Older versions are no longer supported. ** Ticket #199: Particle lifetime cutoff Users can set a maximum lifetime. Particles that live longer than this are not decayed. ** Ticket #215: Madgraph Fixed a problem in reading certain Madgraph event files. ** Les Houches files Optional rescaling of energy or mass of particles which were read in. Check the file to trap 'nan' and 'inf' values early during the read-in. ** File cleanup Most inline functions are now defined in the headers rather than a separate .icc file. * ThePEG-1.3.0 release: 2008-06-20 ** Statistical errors Error estimates on the cross-sections are now reported in the .out files. ** Decaymode setup The decaymode setup has been reworked, keeping backwards compatibility. The 'defaultparticle' command has been removed. ** Madgraph reader Updated to latest Madevent version. ** LHAPDF improvements lhapdf-config is not used anymore to determine the location of the PDF sets. Instead, they are fixed at configure-time. ** HepMC The beam particles are now set correctly. ** Arbitrary search paths for .in files Input files don't have to be in '.', use '-I' to specify additional directories. ** Ticket #172 Fix for cuts in mirrored processes. ** Simpler emacs macros The emacs macros have been significantly cleaned up. ** Helicity vertex classes Potential prolems with uninitialized variables were fixed. No actual bugs had occurred from here. ** Memory leak fixes and performance Several loops of shared pointers were fixed involving DecayModes and ParticleData. FixedSizeAllocator was removed, the regular 'new' and 'delete' is now used for all allocations, giving a 5% speedup. * ThePEG-1.2.0 release: 2008-04-18 ** ThePEG uses GSL The GNU Scientific Library and its headers are now required for building ThePEG. RandomGenerator partially uses GSL now. ** Ticket #160: HepMC converter Optionally, the HepMC converter can fill an external GenEvent object, instead of newing one internally. ** 'globallibrary' command The 'globallibrary' command can be used to register libraries which are useful for the whole run. They do not need to be listed in each class's library() function. ** Resume from dump files (also solves #149) The --resume command line flag instructs runThePEG to resume the run from a previous dump file. 'set Generator:DumpPeriod 1000' writes such a checkpoint every 1000 events. ** New Repository interface (also solves #141) The repository now provides a Repository::cleanup() method, to be called at the end of a run. Alternatively, a Repository object can be created and its member functions called. ** Les Houches interface improvements. LesHouchesReader had had a major overhaul. ** XSecCheck analysis Issues a warning if a target cross section is not met in the run. ** Weighted events Handling of weighted events was improved. ** Ticket #124: 'Deleted' interface option Interfaces can be declared obsolete, a warning will be issued when they're used. ** LHAPDF interface The interface now allows photons as partons inside a hadron. ** gcc 4.3.0 ThePEG compiles cleanly with the new gcc 4.3 series. ** Bug #138 Lepton NoPDF now works. * ThePEG-1.1.2 release: 2008-02-25 ** Bug #136: dSigHatDR Efficiency improvement for zero PDF. ** Bug #137: HepMC conversion The HepMC converter now takes a unit argument to specify which energy and length units should be used in the HepMC event. ** Bug #140 / #141: Crash on shutdown This fixes a bug introduced in 1.1.1. External code interfaces should now work again. ** Bug #151: Loop in Remnant record The loop between remnants in the event record has been removed. ** PDF improvements Fix to handling of maximum flavour from LHAPDF. Fix for xmin calculation. * ThePEG-1.1.1 release: 2007-12-07 ** Bug #46: Reproducibility Fixed a problem where runs were not identical for a given random number seed. You now _must_ reset the seed if you need independent event generator runs. ** Detection of gcc abs bug http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34130 configure now checks for and works around the gcc abs() bug. ** Separate LWHFactory library Fixed problem in Rivet interaction by factoring out LWHFactory into its own dynamic library. * ThePEG-1.1.0 release: 2007-11-20 ** New vector classes ThePEG now uses its own internal Vector classes, making it independent of CLHEP. ** New dimension checking mechanism Optionally, any physical expression is checked for dimensional correctness at compile time. ** Extended Helicity classes A full set of helicity amplitude classes has been transferred from Herwig++. ** unlisted Many other improvements and small bug fixes, see ChangeLog. * ThePEG-1.0.1 release: 2006-11-22 ** unlisted: Fixed memory leak in LesHouchesReader. ** Bug #58 maximumCMEnergy() member of the EventGenerator returns zero. See ChangeLog entry 2006-10-06 ** Bug #62 fixed 'const' behaviour in Lorentz spinor classes ** Bug #68 Improved error message for switch options ** unlisted Improved compile-time LHAPDF library and include file handling. ** unlisted Bug in IteratorRange::rrange(const Container &). ** unlisted fixed Selector::swap() ** unlisted Bug in ClusterCollapser where no colour-singlet particles were considered for momentum compensation if no coloured particles were present. ** unlisted Bug in LeptonLeptonRemnant: minX variable not persistent ** unlisted scale of the produced coloured particles was not set in Onium3GDecayer and ColourPairDecayer ** unlisted unused default path removed from DynamicLoader * ThePEG-1.0 release: 2006-09-27 diff --git a/PDF/LHAPDF6.cc b/PDF/LHAPDF6.cc --- a/PDF/LHAPDF6.cc +++ b/PDF/LHAPDF6.cc @@ -1,478 +1,512 @@ // -*- C++ -*- // // LHAPDF6.cc is a part of ThePEG - Toolkit for HEP Event Generation // // Copyright (C) 2014-2019 Leif Lonnblad, David Grellscheid // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the LHAPDF class. // #include "LHAPDF6.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/Interface/Deleted.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "LHAPDF/LHAPDF.h" using std::vector; using std::string; using std::pair; using ThePEG::GeV2; using ThePEG::cPDVector; ThePEG::LHAPDF::LHAPDF() : thePDF(), thePDFName("THEPEG_NO_PDFSET_CHOSEN"), theMember(0), theMaxFlav(5), xMin(0.0), xMax(1.0), Q2Min(ZERO), Q2Max(Constants::MaxEnergy2) {} ThePEG::LHAPDF::LHAPDF(const LHAPDF & x) : PDFBase(x), thePDF(), thePDFName(x.thePDFName), theMember(x.theMember), theMaxFlav(x.theMaxFlav), xMin(x.xMin), xMax(x.xMax), Q2Min(x.Q2Min), Q2Max(x.Q2Max) {} ThePEG::IBPtr ThePEG::LHAPDF::clone() const { return new_ptr(*this); } ThePEG::IBPtr ThePEG::LHAPDF::fullclone() const { return new_ptr(*this); } void ThePEG::LHAPDF::initPDFptr() { ::LHAPDF::setVerbosity(std::max(0, ThePEG::Debug::level - 1)); if ( thePDF && thePDF->set().name() == thePDFName && thePDF->memberID() == theMember ) return; delete thePDF; thePDF = ::LHAPDF::mkPDF(thePDFName, theMember); xMin = thePDF->xMin(); xMax = thePDF->xMax(); Q2Min = thePDF->q2Min() * GeV2; Q2Max = thePDF->q2Max() * GeV2; } void ThePEG::LHAPDF::doinit() { PDFBase::doinit(); initPDFptr(); } void ThePEG::LHAPDF::dofinish() { PDFBase::dofinish(); delete thePDF; thePDF = 0; } void ThePEG::LHAPDF::doinitrun() { PDFBase::doinitrun(); initPDFptr(); } void ThePEG::LHAPDF::setPDFName(string name) { if ( ::LHAPDF::endswith(name, ".LHgrid") ) { name = name.substr(0, name.size() - 7); } else if ( ::LHAPDF::endswith(name, ".LHpdf") ) { name = name.substr(0, name.size() - 6); } // fix the eternal typo if ( name == "cteq6ll" ) name = "cteq6l1"; if ( ::LHAPDF::contains(::LHAPDF::availablePDFSets(), name) ) { thePDFName = name; theMember = 0; } else { Throw() << "'set " << fullName() << ":PDFName " << name << "': PDF not installed. Try 'lhapdf install'.\n" << Exception::setuperror; } } void ThePEG::LHAPDF::setPDFMember(int member) { try { ::LHAPDF::PDFInfo * test = ::LHAPDF::mkPDFInfo(thePDFName, member); if ( test ) theMember = member; delete test; } catch (::LHAPDF::ReadError & e) { Throw() << e.what() << Exception::setuperror; } } string ThePEG::LHAPDF::doTest(string input) { double x = 0; Energy2 Q2 = ZERO; Energy2 P2 = ZERO; istringstream is(input); is >> x >> iunit(Q2, GeV2) >> iunit(P2, GeV2); initPDFptr(); ostringstream os; for ( int i = 0; i < 13; ++i ) os << " " << thePDF->xfxQ2(i,x,Q2/GeV2); return os.str(); } bool ThePEG::LHAPDF::canHandleParticle(tcPDPtr particle) const { using namespace ParticleID; return abs(particle->id()) == pplus || abs(particle->id()) == n0; } cPDVector ThePEG::LHAPDF::partons(tcPDPtr particle) const { // We assume that all standard partons can be extracted. const ::LHAPDF::PDFSet & pdfset = ::LHAPDF::getPDFSet(thePDFName); const vector & flavs = pdfset.get_entry_as< vector >("Flavors"); - cPDVector ret; - ret.reserve( flavs.size() ); - if ( canHandleParticle(particle) ) { - for ( size_t i=0; i < flavs.size(); ++i ) { - ret.push_back(getParticleData(flavs[i])); + // check for leptons and antileptons + bool hasLepton [3] = {false,false,false}; + bool hasAntiLepton[3] = {false,false,false}; + vector flav2; + flav2.reserve(flavs.size()+3); + for (int val : flavs) { + if(abs(val) == 11 || abs(val) ==13 || abs(val)==15) { + int ilep = (abs(val)-11)/2; + if(val>0) hasLepton[ilep]=true; + else hasAntiLepton[ilep]=true; + } + flav2.push_back(val); + } + for(unsigned int ix=0;ix<3;++ix) { + if(hasLepton[ix] && !hasAntiLepton[ix]) { + int pid = 11+2*ix; + cerr << "Charged lepton with pid=" << pid << " found in PDF but not the antilepton" + << " assumming antilepton PDF is equal to the lepton one\n"; + flav2.push_back(-pid); } } - assert( ret.size() == flavs.size() ); + + cPDVector ret; + ret.reserve( flav2.size() ); + if ( canHandleParticle(particle) ) { + for ( size_t i=0; i < flav2.size(); ++i ) { + ret.push_back(getParticleData(flav2[i])); + } + } + assert( ret.size() == flav2.size() ); return ret; } double ThePEG::LHAPDF::xfx(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale, double x, double, Energy2) const { // Here we should return the actual density. using namespace ThePEG::ParticleID; double Q2 = partonScale/GeV2; if ( ! thePDF->inRangeXQ2(x, Q2) ) { switch ( rangeException ) { case rangeThrow: Throw() << "Momentum fraction (x=" << x << ") or scale (Q2=" << Q2 << " GeV^2) was outside of limits in PDF " << name() << "." << Exception::eventerror; break; case rangeZero: return 0.0; case rangeFreeze: x = min(max(x, xMin), xMax); Q2 = min(max(Q2, Q2Min/GeV2), Q2Max/GeV2); } } int pid = parton->id(); int abspid = abs(pid); switch ( pid ) { case t: case tbar: case b: case bbar: case c: case cbar: return maxFlav() < abspid ? 0.0 : thePDF->xfxQ2(pid,x,Q2); case s: case sbar: return thePDF->xfxQ2(pid,x,Q2); case u: switch ( particle->id() ) { case n0: return thePDF->xfxQ2(d ,x,Q2); case pbarminus: return thePDF->xfxQ2(ubar,x,Q2); case nbar0: return thePDF->xfxQ2(dbar,x,Q2); case pplus: default: return thePDF->xfxQ2(u ,x,Q2); } case ubar: switch ( particle->id() ) { case n0: return thePDF->xfxQ2(dbar,x,Q2); case pbarminus: return thePDF->xfxQ2(u ,x,Q2); case nbar0: return thePDF->xfxQ2(d ,x,Q2); case pplus: default: return thePDF->xfxQ2(ubar,x,Q2); } case d: switch ( particle->id() ) { case n0: return thePDF->xfxQ2(u ,x,Q2); case pbarminus: return thePDF->xfxQ2(dbar,x,Q2); case nbar0: return thePDF->xfxQ2(ubar,x,Q2); case pplus: default: return thePDF->xfxQ2(d ,x,Q2); } case dbar: switch ( particle->id() ) { case n0: return thePDF->xfxQ2(ubar,x,Q2); case pbarminus: return thePDF->xfxQ2(d ,x,Q2); case nbar0: return thePDF->xfxQ2(u ,x,Q2); case pplus: default: return thePDF->xfxQ2(dbar,x,Q2); } case g: return thePDF->xfxQ2(g,x,Q2); case ParticleID::gamma: return thePDF->xfxQ2(ParticleID::gamma,x,Q2); + case ParticleID::eminus: case ParticleID::muminus: case ParticleID::tauminus: + case ParticleID::eplus : case ParticleID::muplus : case ParticleID::tauplus : + if(pid>0 || thePDF->hasFlavor(pid)) + return thePDF->xfxQ2(pid,x,Q2); + else + return thePDF->xfxQ2(-pid,x,Q2); } return 0.0; } double ThePEG::LHAPDF::xfvl(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale, double l, Energy2 particleScale) const { using Math::exp1m; return xfvx(particle, parton, partonScale, exp(-l), exp1m(-l), particleScale); } double ThePEG::LHAPDF::xfvx(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale, double x, double, Energy2) const { // Here we should return the actual valence density. This will only // work properly for nucleons using namespace ThePEG::ParticleID; double Q2 = partonScale / GeV2; if ( ! thePDF->inRangeXQ2(x, Q2) ) { switch ( rangeException ) { case rangeThrow: Throw() << "Momentum fraction (x=" << x << ") or scale (Q2=" << Q2 << " GeV^2) was outside of limits in PDF " << name() << "." << Exception::eventerror; break; case rangeZero: return 0.0; case rangeFreeze: x = min(max(x, xMin), xMax); Q2 = min(max(Q2, Q2Min/GeV2), Q2Max/GeV2); } } switch ( parton->id() ) { case t: case tbar: case b: case bbar: case c: case cbar: case s: case sbar: case ParticleID::gamma: return 0.0; case u: switch ( particle->id() ) { case n0: return thePDF->xfxQ2(d,x,Q2) - thePDF->xfxQ2(dbar,x,Q2); case pbarminus: return 0.0; case nbar0: return 0.0; case pplus: return thePDF->xfxQ2(u,x,Q2) - thePDF->xfxQ2(ubar,x,Q2); default: return 0.0; } case ubar: switch ( particle->id() ) { case n0: return 0.0; case pbarminus: return thePDF->xfxQ2(u,x,Q2) - thePDF->xfxQ2(ubar,x,Q2); case nbar0: return thePDF->xfxQ2(d,x,Q2) - thePDF->xfxQ2(dbar,x,Q2); case pplus: default: return 0.0; } case d: switch ( particle->id() ) { case n0: return thePDF->xfxQ2(u,x,Q2) - thePDF->xfxQ2(ubar,x,Q2); case pbarminus: return 0.0; case nbar0: return 0.0; case pplus: return thePDF->xfxQ2(d,x,Q2) - thePDF->xfxQ2(dbar,x,Q2); default: return 0.0; } case dbar: switch ( particle->id() ) { case n0: return 0.0; case pbarminus: return thePDF->xfxQ2(d,x,Q2) - thePDF->xfxQ2(dbar,x,Q2); case nbar0: return thePDF->xfxQ2(u,x,Q2) - thePDF->xfxQ2(ubar,x,Q2); case pplus: default: return 0.0; } case ParticleID::g: return 0.0; } return 0.0; } double ThePEG::LHAPDF::xfsx(tcPDPtr particle, tcPDPtr parton, Energy2 partonScale, double x, double, Energy2) const { // Here we should return the actual density. using namespace ThePEG::ParticleID; double Q2 = partonScale / GeV2; if ( ! thePDF->inRangeXQ2(x, Q2) ) { switch ( rangeException ) { case rangeThrow: Throw() << "Momentum fraction (x=" << x << ") or scale (Q2=" << Q2 << " GeV^2) was outside of limits in PDF " << name() << "." << Exception::eventerror; break; case rangeZero: return 0.0; case rangeFreeze: x = min(max(x, xMin), xMax); Q2 = min(max(Q2, Q2Min/GeV2), Q2Max/GeV2); } } int pid = parton->id(); int abspid = abs(pid); switch ( pid ) { case t: case tbar: case b: case bbar: case c: case cbar: return maxFlav() < abspid ? 0.0 : thePDF->xfxQ2(pid,x,Q2); case s: case sbar: return thePDF->xfxQ2(pid,x,Q2); case u: switch ( particle->id() ) { case n0: return thePDF->xfxQ2(dbar,x,Q2); case pbarminus: return thePDF->xfxQ2(ubar,x,Q2); case nbar0: return thePDF->xfxQ2(dbar,x,Q2); case pplus: return thePDF->xfxQ2(ubar,x,Q2); default: return thePDF->xfxQ2(u ,x,Q2); } case ubar: switch ( particle->id() ) { case n0: return thePDF->xfxQ2(dbar,x,Q2); case pbarminus: return thePDF->xfxQ2(ubar,x,Q2); case nbar0: return thePDF->xfxQ2(dbar,x,Q2); case pplus: default: return thePDF->xfxQ2(ubar,x,Q2); } case d: switch ( particle->id() ) { case n0: return thePDF->xfxQ2(ubar,x,Q2); case pbarminus: return thePDF->xfxQ2(dbar,x,Q2); case nbar0: return thePDF->xfxQ2(ubar,x,Q2); case pplus: return thePDF->xfxQ2(dbar,x,Q2); default: return thePDF->xfxQ2(d ,x,Q2); } case dbar: switch ( particle->id() ) { case n0: return thePDF->xfxQ2(ubar,x,Q2); case pbarminus: return thePDF->xfxQ2(dbar,x,Q2); case nbar0: return thePDF->xfxQ2(ubar,x,Q2); case pplus: default: return thePDF->xfxQ2(dbar,x,Q2); } case g: return thePDF->xfxQ2(g,x,Q2); case ParticleID::gamma: return thePDF->xfxQ2(ParticleID::gamma,x,Q2); + case ParticleID::eminus: case ParticleID::muminus: case ParticleID::tauminus: + case ParticleID::eplus : case ParticleID::muplus : case ParticleID::tauplus : + if(pid>0 || thePDF->hasFlavor(pid)) + return thePDF->xfxQ2(pid,x,Q2); + else + return thePDF->xfxQ2(-pid,x,Q2); } return 0.0; } void ThePEG::LHAPDF::persistentOutput(PersistentOStream & os) const { os << thePDFName << theMember << theMaxFlav << xMin << xMax << ounit(Q2Min, GeV2) << ounit(Q2Max, GeV2); } void ThePEG::LHAPDF::persistentInput(PersistentIStream & is, int) { is >> thePDFName >> theMember >> theMaxFlav >> xMin >> xMax >> iunit(Q2Min, GeV2) >> iunit(Q2Max, GeV2); initPDFptr(); } ThePEG::DescribeClass describeLHAPDF("ThePEG::LHAPDF", "ThePEGLHAPDF.so"); void ThePEG::LHAPDF::Init() { static ClassDocumentation documentation ("The LHAPDF class inherits from PDFBase and implements an interface " "to the LHAPDF library of parton density function parameterizations. " "This class is available even if LHAPDF was not properly installed " "when ThePEG was installed, but will then produce an error in the " "initialization. Note that the valence densities from the xfvx() and " "xfvl() function will only work properly for nucleons. All other " "particles will have zero valence densities."); static Deleted interfacePType ("PType", "The LHAPDFv6 interface currently does not support pi."); static Parameter interfacePDFName ("PDFName", "The name of the PDF set to be used. Should correspond to " "the LHAPDF v6 name.", &ThePEG::LHAPDF::thePDFName, "THEPEG_NO_PDFSET_CHOSEN", true, false, &ThePEG::LHAPDF::setPDFName); static Deleted interfacePDFNumber ("PDFNumber", "Not implemented in the LHAPDFv6 interface. " "Use :PDFName and :Member instead."); static Parameter interfaceMember ("Member", "The chosen member of the selected PDF set.", &ThePEG::LHAPDF::theMember, 0, 0, 0, true, false, Interface::lowerlim, &ThePEG::LHAPDF::setPDFMember); static Deleted interfacePDFLIBNumbers ("PDFLIBNumbers", "Not implemented in the LHAPDFv6 interface. " "Use :PDFName and :Member instead."); static Deleted interfaceEnablePartonicGamma ("EnablePartonicGamma", "Not required in LHAPDFv6."); static Deleted interfacePhotonOption ("PhotonOption", "Not required in LHAPDFv6."); static Deleted interfaceMaxNSet ("MaxNSet", "Not required in LHAPDFv6."); static Command interfaceTest ("Test", "Write out the values of the chosen PDF set using the x, Q2 and P2 " "parameters supplied.", &ThePEG::LHAPDF::doTest, true); static Deleted interfaceVerboseLevel ("VerboseLevel", "LHAPDFv6 uses general debug level instead."); static Parameter interfaceMaxFlav ("MaxFlav", "The maximum number of flavours for which non-zero densities are " "reported. The actual number of flavours may be less depending on " "the chosen PDF set.", &ThePEG::LHAPDF::theMaxFlav, 5, 3, 0, true, false, Interface::lowerlim); interfacePDFName.rank(9); interfaceMember.rank(8); }