diff --git a/MatrixElement/FxFx/FxFxFileReader.cc b/MatrixElement/FxFx/FxFxFileReader.cc --- a/MatrixElement/FxFx/FxFxFileReader.cc +++ b/MatrixElement/FxFx/FxFxFileReader.cc @@ -1,914 +1,929 @@ // -*- C++ -*- // // FxFxFileReader.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 FxFxFileReader class. // #include "FxFxFileReader.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; FxFxFileReader:: FxFxFileReader(const FxFxFileReader & x) : FxFxReader(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) {} FxFxFileReader::~FxFxFileReader() {} IBPtr FxFxFileReader::clone() const { return new_ptr(*this); } IBPtr FxFxFileReader::fullclone() const { return new_ptr(*this); } bool FxFxFileReader::preInitialize() const { return true; } void FxFxFileReader::doinit() { FxFxReader::doinit(); // 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() << "FxFxFileReader::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() << "FxFxFileReader::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 FxFxFileReader::doinit()." << " Please check the Les Houches file.\n" << Exception::runerror; tcPDPtr p = getParticleData(t); if( !p ) throw SetupException() << "FxFxFileReader::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() << "FxFxFileReader::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 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() << "FxFxFileReader::doinit() Decayer must be set using the " << "FxFxFileReader:Decayer" << " must be set to allow the creation of new" << " decay modes." << Exception::runerror; if(!dm) { dm = generator()->preinitCreateDecayMode(tag); if(!dm) Throw() << "FxFxFileReader::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 FxFxFileReader::optWeightNamesFunc() { return optionalWeightsNames; } vector FxFxFileReader::optWeightsNamesFunc() { return optionalWeightsNames; } void FxFxFileReader::open() { if ( filename().empty() ) throw FxFxFileError() << "No Les Houches file name. " << "Use 'set " << name() << ":FileName'." << Exception::runerror; cfile.open(filename()); if ( !cfile ) throw FxFxFileError() << "The FxFxFileReader '" << 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"]; //cout << LHFVersion << endl; 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; } ++ws; } while (isc); /* now get the relevant information * e.g. scales or PDF sets used */ string startDEL = ">"; //starting delimiter string stopDEL = ""; //end delimiter unsigned firstLim = hs.find(startDEL); //find start of delimiter // 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); + //cout << "scinfo = " << scinfo << endl; /* 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; string ecomstring = "ecom"; optionalWeights[ecomstring.c_str()] = heprup.EBMUP.first+heprup.EBMUP.second; 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; 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 + //cout << "weightName, weightValue= " << weightName << ", " << weightValue << endl; 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); string str_newline= "\n"; erase_substr(mg5clusinfo,str_newline); 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){ + string itfirst = it->first; + erase_substr(itfirst, "'"); + erase_substr(itfirst, "\""); for (map::const_iterator it2=scalemap.begin(); it2!=scalemap.end(); ++it2){ //find the scale id in the scale information and add this information - if(it->first==it2->first) { + string it2first = it2->first; + erase_substr(it2first, "'"); + erase_substr(it2first, "\""); + //cout << "itfirst, it2first = " << itfirst << "\t" << it2first << endl; + if(itfirst==it2first) { 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 FxFxFileReader::close() { cfile.close(); } void FxFxFileReader::persistentOutput(PersistentOStream & os) const { os << neve << LHFVersion << outsideBlock << headerBlock << initComments << initAttributes << eventComments << eventAttributes << theFileName << theQNumbers << theIncludeFxFxTags << theIncludeCentral << theDecayer; } void FxFxFileReader::persistentInput(PersistentIStream & is, int) { is >> neve >> LHFVersion >> outsideBlock >> headerBlock >> initComments >> initAttributes >> eventComments >> eventAttributes >> theFileName >> theQNumbers >> theIncludeFxFxTags >> theIncludeCentral >> theDecayer; ieve = 0; } ClassDescription FxFxFileReader::initFxFxFileReader; // Definition of the static class description member. void FxFxFileReader::Init() { static ClassDocumentation documentation ("ThePEG::FxFxFileReader 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.", &FxFxFileReader::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.", &FxFxFileReader::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", &FxFxFileReader::theIncludeFxFxTags, true, 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", &FxFxFileReader::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", &FxFxFileReader::theDecayer, false, false, true, true, false); } void FxFxFileReader::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/Models/Feynrules/python/ufo2herwig b/Models/Feynrules/python/ufo2herwig old mode 100755 new mode 100644 --- a/Models/Feynrules/python/ufo2herwig +++ b/Models/Feynrules/python/ufo2herwig @@ -1,474 +1,475 @@ #! /usr/bin/env python from __future__ import division from __future__ import print_function import os, sys, pprint, argparse, re,copy from string import Template # add path to the ufo conversion modules modulepath = os.path.join("@PKGLIBDIR@",'python') sys.path.append(modulepath) from ufo2peg import * # set up the option parser for command line input parser = argparse.ArgumentParser( description='Create Herwig model files from Feynrules UFO input.' ) parser.add_argument( 'ufodir', metavar='UFO_directory', help='the UFO model directory' ) parser.add_argument( '-v', '--verbose', action="store_true", help="print verbose output" ) parser.add_argument( '-n','--name', default="FRModel", help="set custom nametag for the model" ) parser.add_argument( '--ignore-skipped', action="store_true", help="ignore skipped vertices and produce output anyway" ) parser.add_argument( '--split-model', action="store_true", help="Split the model file into pieces to improve compilation for models with many parameters" ) parser.add_argument( '--no-generic-loop-vertices', action="store_true", help="Don't include the automatically generated generic loop vertices for h->gg and h->gamma gamma" ) parser.add_argument( '--include-generic', action="store_true", help="Include support for generic spin structures" ) parser.add_argument( '--use-Herwig-Higgs', action="store_true", help="Use the internal Herwig SM modeling and vertices for Higgs interactions, there may be sign issues but some UFO models have very minimal Higgs interactions" ) parser.add_argument( '--use-generic-for-tensors', action="store_true", help="Use the generic machinery for all tensor vertices (debugging only)" ) parser.add_argument( '--forbidden-particle-name', action="append", default=["eta","phi"], help="Add particle names not allowed as names in UFO models to avoid conflicts with"+\ "Herwig internal particles, names will have _UFO appended" ) parser.add_argument( '--convert', action="store_true", default=False, help="Convert the UFO model for python2 to python3 before loading it." ) # get the arguments args = parser.parse_args() # convert model to python 3 if needed if(args.convert) : + prepForConversion(args.ufodir) convertToPython3(args.ufodir) # import the model try : import importlib,importlib.util try : path,mod = os.path.split(os.path.abspath(args.ufodir)) spec = importlib.util.spec_from_file_location(mod,os.path.join(os.path.abspath(args.ufodir),"__init__.py")) FR = spec.loader.load_module() except : print ("Could not load the UFO python module.\n", "This is usually because you are using python3 and the UFO model is in python2.\n", "The --convert option can be used to try and convert it, but there is no guarantee.\n", "As this modifies the model you should make sure you have a backup copy.\n") quit() except : print ("Newer python import did not work, reverting to python 2 imp module") import imp path,mod = os.path.split(os.path.abspath(args.ufodir)) FR = imp.load_module(mod,*imp.find_module(mod,[path])) ################################################## ################################################## # get the Model name from the arguments modelname = args.name libname = modelname + '.so' # define arrays and variables #allplist = "" parmdecls = [] parmgetters = [] parmconstr = [] doinit = [] paramstoreplace_ = [] paramstoreplace_expressions_ = [] # get external parameters for printing parmsubs = dict( [ (p.name, p.value) for p in FR.all_parameters if p.nature == 'external' ] ) # evaluate python cmath def evaluate(x): import cmath return eval(x, {'cmath':cmath, 'complexconjugate':FR.function_library.complexconjugate, 'im':FR.function_library.im, 're':FR.function_library.re}, parmsubs) ## get internal params into arrays internal = ( p for p in FR.all_parameters if p.nature == 'internal' ) #paramstoreplaceEW_ = [] #paramstoreplaceEW_expressions_ = [] # calculate internal parameters for p in internal: parmsubs.update( { p.name : evaluate(p.value) } ) # if 'aS' in p.value and p.name != 'aS': # paramstoreplace_.append(p.name) # paramstoreplace_expressions_.append(p.value) # if 'aEWM1' in p.value and p.name != 'aEWM1': # paramstoreplaceEW_.append(p.name) # paramstoreplaceEW_expressions_.append(p.value) parmvalues=copy.copy(parmsubs) # more arrays used for substitution in templates paramsforstream = [] parmmodelconstr = [] # loop over parameters and fill in template stuff according to internal/external and complex/real # WARNING: Complex external parameter input not tested! if args.verbose: print ('verbose mode on: printing all parameters') print ('-'*60) paramsstuff = ('name', 'expression', 'default value', 'nature') pprint.pprint(paramsstuff) interfacedecl_T = """\ static Parameter<{modelname}, {type}> interface{pname} ("{pname}", "The interface for parameter {pname}", &{modelname}::{pname}_, {value}, 0, 0, false, false, Interface::nolimits); """ # sort out the couplings couplingDefns = { "QED" : 99, "QCD" : 99 } try : for coupling in FR.all_orders: name = coupling.name.upper() couplingDefns[name]= coupling.expansion_order except: for coupling in FR.all_couplings: for name,value in coupling.order.items(): if(name not in couplingDefns) : couplingDefns[name]=99 # sort out the particles massnames = {} widthnames = {} for particle in FR.all_particles: # skip ghosts and goldstones if(isGhost(particle) or isGoldstone(particle)) : continue if particle.mass != 'ZERO' and particle.mass.name != 'ZERO': if(particle.mass in massnames) : if(abs(particle.pdg_code) not in massnames[particle.mass]) : massnames[particle.mass].append(abs(particle.pdg_code)) else : massnames[particle.mass] = [abs(particle.pdg_code)] if particle.width != 'ZERO' and particle.width.name != 'ZERO': if(particle.width in widthnames) : if(abs(particle.pdg_code) not in widthnames[particle.width]) : widthnames[particle.width].append(abs(particle.pdg_code)) else : widthnames[particle.width] = [abs(particle.pdg_code)] interfaceDecls = [] modelparameters = {} for p in FR.all_parameters: value = parmsubs[p.name] if p.type == 'real': assert( value.imag < 1.0e-16 ) value = value.real if p.nature == 'external': if p not in massnames and p not in widthnames: interfaceDecls.append( interfacedecl_T.format(modelname=modelname, pname=p.name, value=value, type=typemap(p.type)) ) else: interfaceDecls.append('\n// no interface for %s. Use particle definition instead.\n' % p.name) if hasattr(p,'lhablock'): lhalabel = '{lhablock}_{lhacode}'.format( lhablock=p.lhablock.upper(), lhacode='_'.join(map(str,p.lhacode)) ) if p not in massnames and p not in widthnames: parmmodelconstr.append('set %s:%s ${%s}' % (modelname, p.name, lhalabel)) else: parmmodelconstr.append('# %s is taken from the particle setup' % p.name) modelparameters[lhalabel] = value parmsubs[p.name] = lhalabel else: if p not in massnames and p not in widthnames: parmmodelconstr.append('set %s:%s %s' % (modelname, p.name, value)) else: parmmodelconstr.append('# %s is taken from the particle setup' % p.name) parmsubs[p.name] = value if p not in massnames and p not in widthnames: parmconstr.append('%s_(%s)' % (p.name, value)) else: parmconstr.append('%s_()' % p.name) else : parmconstr.append('%s_()' % p.name) parmsubs[p.name] = value elif p.type == 'complex': value = complex(value) if p.nature == 'external': parmconstr.append('%s_(%s,%s)' % (p.name, value.real, value.imag)) else : parmconstr.append('%s_(%s,%s)' % (p.name, 0.,0.)) parmsubs[p.name] = value else: raise Exception('Unknown data type "%s".' % p.type) parmdecls.append(' %s %s_;' % (typemap(p.type), p.name)) parmgetters.append(' %s %s() const { return %s_; }' % (typemap(p.type),p.name, p.name)) paramsforstream.append('%s_' % p.name) expression, symbols = 'return %s_' % p.name, None if p.nature != 'external': expression, symbols = py2cpp(p.value) text = add_brackets(expression, symbols) text=text.replace('()()','()') doinit.append(' %s_ = %s;' % (p.name, text) ) if(p.type == 'complex') : doinit.append(' if(std::isnan(%s_.real()) || std::isnan(%s_.imag()) || std::isinf(%s_.real()) || std::isinf(%s_.imag())) {throw InitException() << "Calculated parameter %s is nan or inf in Feynrules model. Check your input parameters.";}' % (p.name,p.name,p.name,p.name,p.name) ) else : doinit.append(' if(std::isnan(%s_) || std::isinf(%s_)) {throw InitException() << "Calculated parameter %s is nan or inf in Feynrules model. Check your input parameters.";}' % (p.name,p.name,p.name) ) if p in massnames: for idCode in massnames[p] : doinit.append(' resetMass(%s,%s_ * GeV);' % (idCode, p.name) ) if p in widthnames: for idCode in widthnames[p] : doinit.append(' getParticleData(%s)->width(%s_ * GeV);' % (idCode, p.name) ) doinit.append(' getParticleData(%s)->cTau (%s_ == 0.0 ? Length() : hbarc/(%s_*GeV));' % (idCode, p.name, p.name) ) doinit.append(' getParticleData(%s)->widthCut(10. * %s_ * GeV);' % (idCode, p.name) ) elif p.nature == 'external': if p in massnames: for idCode in massnames[p] : doinit.append(' %s_ = getParticleData(%s)->mass() / GeV;' % (p.name, idCode) ) if p in widthnames: for idCode in widthnames[p] : doinit.append(' %s_ = getParticleData(%s)->width() / GeV;' % (p.name, idCode) ) if args.verbose: pprint.pprint((p.name,p.value, value, p.nature)) pcwriter = ParamCardWriter(FR.all_parameters) paramcard_output = '\n'.join(pcwriter.output) ### special treatment # if p.name == 'aS': # expression = '0.25 * sqr(strongCoupling(q2)) / Constants::pi' # elif p.name == 'aEWM1': # expression = '4.0 * Constants::pi / sqr(electroMagneticCoupling(q2))' # elif p.name == 'Gf': # expression = 'generator()->standardModel()->fermiConstant() * GeV2' paramconstructor=': ' for ncount in range(0,len(parmconstr)) : paramconstructor += parmconstr[ncount] if(ncount != len(parmconstr) -1) : paramconstructor += ',' if(ncount%5 == 0 ) : paramconstructor += "\n" paramout="" paramin ="" for ncount in range(0,len(paramsforstream)) : if(ncount !=0 ) : paramout += "<< " + paramsforstream[ncount] paramin += ">> " + paramsforstream[ncount] else : paramout += paramsforstream[ncount] paramin += paramsforstream[ncount] if(ncount%5 == 0 ) : paramout += "\n" paramin += "\n" parmtextsubs = { 'parmgetters' : '\n'.join(parmgetters), 'parmdecls' : '\n'.join(parmdecls), 'parmconstr' : paramconstructor, 'getters' : '', 'decls' : '', 'addVertex' : '', 'doinit' : '\n'.join(doinit), 'ostream' : paramout, 'istream' : paramin , 'refs' : '', 'parmextinter': ''.join(interfaceDecls), 'ModelName': modelname, 'calcfunctions': '', 'param_card_data': paramcard_output } ################################################## ################################################## ################################################## # set up the conversion of the vertices vertexConverter = VertexConverter(FR,parmvalues,couplingDefns) vertexConverter.readArgs(args) # convert the vertices vertexConverter.convert() cdefs="" couplingOrders="" ncount=2 for name,value in couplingDefns.items() : if(name=="QED") : couplingOrders+=" setCouplings(\"%s\",make_pair(%s,%s));\n" %(name,1,value) elif (name=="QCD") : couplingOrders+=" setCouplings(\"%s\",make_pair(%s,%s));\n" %(name,2,value) else : ncount+=1 cdefs +=" const T %s = %s;\n" % (name,ncount) couplingOrders+=" setCouplings(\"%s\",make_pair(%s,%s));\n" % (name,ncount,value) # coupling definitions couplingTemplate= """\ namespace ThePEG {{ namespace Helicity {{ namespace CouplingType {{ typedef unsigned T; /** * Enums for couplings */ {coup} }} }} }} """ if(cdefs!="") : cdefs = couplingTemplate.format(coup=cdefs) parmtextsubs['couplings'] = cdefs parmtextsubs['couplingOrders'] = couplingOrders # particles plist, names = thepeg_particles(FR,parmsubs,modelname,modelparameters,args.forbidden_particle_name,args.use_Herwig_Higgs) particlelist = [ "# insert HPConstructor:Outgoing 0 /Herwig/{n}/Particles/{p}".format(n=modelname,p=p) for p in names ] # make the first one active to have something runnable in the example .in file particlelist[0] = particlelist[0][2:] particlelist = '\n'.join(particlelist) modelfilesubs = { 'plist' : plist, 'vlist' : vertexConverter.get_vertices(libname), 'setcouplings': '\n'.join(parmmodelconstr), 'ModelName': modelname } # write the files from templates according to the above subs if vertexConverter.should_print(): MODEL_HWIN = getTemplate('LHC-FR.in') if(not args.split_model) : MODEL_CC = [getTemplate('Model.cc')] else : MODEL_EXTRA_CC=getTemplate('Model6.cc') extra_names=[] extra_calcs=[] parmtextsubs['doinit']="" for i in range(0,len(doinit)) : if( i%20 == 0 ) : function_name = "initCalc" +str(int(i/20)) parmtextsubs['doinit'] += function_name +"();\n" parmtextsubs['calcfunctions'] += "void " + function_name + "();\n" extra_names.append(function_name) extra_calcs.append("") extra_calcs[-1] += doinit[i] + "\n" for i in range(0,len(extra_names)) : ccname = '%s_extra_%s.cc' % (modelname,i) writeFile( ccname, MODEL_EXTRA_CC.substitute({'ModelName' : modelname, 'functionName' : extra_names[i], 'functionCalcs' : extra_calcs[i] }) ) MODEL_CC = [getTemplate('Model1.cc'),getTemplate('Model2.cc'),getTemplate('Model3.cc'), getTemplate('Model4.cc'),getTemplate('Model5.cc')] MODEL_H = getTemplate('Model.h') print ('LENGTH',len(MODEL_CC)) MODELINFILE = getTemplate('FR.model') writeFile( 'LHC-%s.in' % modelname, MODEL_HWIN.substitute({ 'ModelName' : modelname, 'Particles' : particlelist }) ) modeltemplate = Template( MODELINFILE.substitute(modelfilesubs) ) writeFile( '%s.h' % modelname, MODEL_H.substitute(parmtextsubs) ) for i in range(0,len(MODEL_CC)) : if(len(MODEL_CC)==1) : ccname = '%s.cc' % modelname else : ccname = '%s.cc' % (modelname + str(i)) writeFile( ccname, MODEL_CC[i].substitute(parmtextsubs) ) writeFile( modelname +'.template', modeltemplate.template ) writeFile( modelname +'.model', modeltemplate.substitute( modelparameters ) ) # copy the Makefile-FR to current directory, # replace with the modelname for compilation with open(os.path.join(modulepath,'Makefile-FR'),'r') as orig: with open('Makefile','w') as dest: dest.write(orig.read().replace("FeynrulesModel.so", libname)) print('finished generating model:\t', modelname) print('model directory:\t\t', args.ufodir) print('generated:\t\t\t', len(FR.all_vertices), 'vertices') print('='*60) print('library:\t\t\t', libname) print('input file:\t\t\t', 'LHC-' + modelname +'.in') print('model file:\t\t\t', modelname +'.model') print('='*60) print("""\ To complete the installation, compile by typing "make". An example input file is provided as LHC-FRModel.in, you'll need to change the required particles in there. """) print('DONE!') print('='*60) diff --git a/Models/Feynrules/python/ufo2peg/__init__.py b/Models/Feynrules/python/ufo2peg/__init__.py --- a/Models/Feynrules/python/ufo2peg/__init__.py +++ b/Models/Feynrules/python/ufo2peg/__init__.py @@ -1,9 +1,10 @@ from .particles import thepeg_particles +from .helpers import prepForConversion from .helpers import convertToPython3 from .helpers import CheckUnique, getTemplate, writeFile, def_from_model from .helpers import add_brackets, typemap,isGhost,isGoldstone from .converter import py2cpp from .write_paramcard import ParamCardWriter from .vertices import VertexConverter diff --git a/Models/Feynrules/python/ufo2peg/helpers.py b/Models/Feynrules/python/ufo2peg/helpers.py --- a/Models/Feynrules/python/ufo2peg/helpers.py +++ b/Models/Feynrules/python/ufo2peg/helpers.py @@ -1,312 +1,338 @@ from __future__ import print_function from string import Template from os import path import sys,cmath,glob import re """ Helper functions for the Herwig Feynrules converter """ class CheckUnique: """Uniqueness checker. - + An object of this class remembers the value it was called with first. Any subsequent call to it will only succeed if the same value is passed again. For example, >>> f = CheckUnique() >>> f(5) >>> f(5) >>> f(4) Traceback (most recent call last): ... AssertionError """ def __init__(self): self.val = None def __call__(self,val): """Store value on first call, then compare.""" if self.val is None: self.val = val else: assert( val == self.val ) def is_number(s): """Check if a value is a number.""" try: float(s) except ValueError: return False return True def getTemplate(name): """Create a python string template from a file.""" templatename = '{name}.template'.format(name=name) # assumes the template files sit next to this script moduledir = path.dirname(path.abspath(__file__)) templatepath = path.join(moduledir,templatename) with open(templatepath, 'r') as f: templateText = f.read() return Template( templateText ) def writeFile(filename, text): """Write text to a filename.""" with open(filename,'w') as f: f.write(text) def coupling_orders(vertex, coupling, defns): # if more than one take QCD over QED and then lowest order in QED if type(coupling) is list: print('not sure this happens') quit() qed = 0 qcd = 0 for coup in coupling : qed1 = coup.order.get('QED',0) qcd1 = coup.order.get('QCD',0) if qcd1 != 0: if qcd == 0 or (qcd1 != 0 and qcd1 < qcd): qcd=qcd1 qed=qed1 else: if qed == 0 or (qed1 != 0 and qed1 < qed): qed=qed1 else: output={} for ctype in defns : output[ctype]=coupling.order.get(ctype,0) return output def def_from_model(FR,s): """Return a C++ line that defines parameter s as coming from the model file.""" if("hw_kine" in s) :return "" stype = typemap(getattr(FR.parameters,s).type) return '{t} {s} = model_->{s}();'.format(t=stype,s=s) _typemap = {'complex':'Complex', 'real':'double'} def typemap(s): return _typemap[s] def add_brackets(expr, syms): result = expr for s in syms: pattern = r'({symb})(\W|$)'.format(symb=s) result = re.sub(pattern, r'\1()\2', result) return result def banner(): return """\ =============================================================================================================== -______ ______ _ __ _ _ _ -| ___| | ___ \ | | / /| | | | (_) _ _ -| |_ ___ _ _ _ __ | |_/ /_ _ | | ___ ___ / / | |_| | ___ _ __ __ __ _ __ _ _| |_ _| |_ +______ ______ _ __ _ _ _ +| ___| | ___ \ | | / /| | | | (_) _ _ +| |_ ___ _ _ _ __ | |_/ /_ _ | | ___ ___ / / | |_| | ___ _ __ __ __ _ __ _ _| |_ _| |_ | _|/ _ \| | | || \_ \ | /| | | || | / _ \/ __| / / | _ | / _ \| \__|\ \ /\ / /| | / _` ||_ _||_ _| -| | | __/| |_| || | | || |\ \| |_| || || __/\__ \ / / | | | || __/| | \ V V / | || (_| | |_| |_| -\_| \___| \__, ||_| |_|\_| \_|\__,_||_| \___||___//_/ \_| |_/ \___||_| \_/\_/ |_| \__, | - __/ | __/ | - |___/ |___/ +| | | __/| |_| || | | || |\ \| |_| || || __/\__ \ / / | | | || __/| | \ V V / | || (_| | |_| |_| +\_| \___| \__, ||_| |_|\_| \_|\__,_||_| \___||___//_/ \_| |_/ \___||_| \_/\_/ |_| \__, | + __/ | __/ | + |___/ |___/ =============================================================================================================== generating model/vertex/.model/.in files please be patient! =============================================================================================================== """ #################### ??? ####################### # function that replaces alphaS (aS)-dependent variables # with their explicit form which also contains strongCoupling def aStoStrongCoup(stringin, paramstoreplace, paramstoreplace_expressions): #print stringin for xx in range(0,len(paramstoreplace)): #print paramstoreplace[xx], paramstoreplace_expressions[xx] stringout = stringin.replace(paramstoreplace[xx], '(' + PyMathToThePEGMath(paramstoreplace_expressions[xx],allparams) + ')') stringout = stringout.replace('aS', '(sqr(strongCoupling(q2))/(4.0*Constants::pi))') #print 'resulting string', stringout return stringout # function that replaces alphaEW (aEW)-dependent variables # with their explicit form which also contains weakCoupling def aEWtoWeakCoup(stringin, paramstoreplace, paramstoreplace_expressions): #print stringin for xx in range(0,len(paramstoreplace)): #print paramstoreplace[xx], paramstoreplace_expressions[xx] stringout = stringin.replace(paramstoreplace[xx], '(' + PyMathToThePEGMath(paramstoreplace_expressions[xx],allparams) + ')') stringout = stringout.replace('aEWM1', '(1/(sqr(electroMagneticCoupling(q2))/(4.0*Constants::pi)))') #print 'resulting string', stringout return stringout if __name__ == "__main__": import doctest doctest.testmod() if False: - + # Check if the Vertex is self-conjugate or not pdgcode = [0,0,0,0] notsmvertex = False vhasw = 0 vhasz = 0 vhasf = 0 vhasg = 0 vhash = 0 vhasp = 0 # print 'printing particles in vertex' for i in range(len(v.particles)): # print v.particles[i].pdg_code pdgcode[i] = v.particles[i].pdg_code if pdgcode[i] == 23: vhasz += 1 if pdgcode[i] == 22: vhasp += 1 if pdgcode[i] == 25: vhash += 1 if pdgcode[i] == 21: vhasg += 1 if pdgcode[i] == 24: vhasw += 1 if abs(pdgcode[i]) < 7 or (abs(pdgcode[i]) > 10 and abs(pdgcode[i]) < 17): vhasf += 1 if pdgcode[i] not in SMPARTICLES: notsmvertex = True - -# treat replacement of SM vertices with BSM vertices? + +# treat replacement of SM vertices with BSM vertices? if notsmvertex == False: if( (vhasf == 2 and vhasz == 1) or (vhasf == 2 and vhasw == 1) or (vhasf == 2 and vhash == 1) or (vhasf == 2 and vhasg == 1) or (vhasf == 2 and vhasp == 0) or (vhasg == 3) or (vhasg == 4) or (vhasw == 2 and vhash == 1) or (vhasw == 3) or (vhasw == 4) or (vhash == 1 and vhasg == 2) or (vhash == 1 and vhasp == 2)): #print 'VERTEX INCLUDED IN STANDARD MODEL!' v.include = 0 notincluded += 1 #continue - - + + selfconjugate = 0 for j in range(len(pdgcode)): for k in range(len(pdgcode)): if j != k and j != 0 and abs(pdgcode[j]) == abs(pdgcode[k]): selfconjugate = 1 #print 'self-conjugate vertex' # print pdgcode[j] # if the Vertex is not self-conjugate, then add the conjugate vertex # automatically scfac = [1,1,1,1] if selfconjugate == 0: #first find the self-conjugate particles for u in range(len(v.particles)): if v.particles[u].selfconjugate == 0: scfac[u] = -1 # print 'particle ', v.particles[u].pdg_code, ' found not to be self-conjugate' - + if selfconjugate == 0: plistarray[1] += str(scfac[1] * v.particles[1].pdg_code) + ',' + str(scfac[0] * v.particles[0].pdg_code) + ',' + str(scfac[2] * v.particles[2].pdg_code) - if len(v.particles) is 4: + if len(v.particles) == 4: plistarray[1] += ',' + str(scfac[3] * v.particles[3].pdg_code) #print 'Conjugate vertex:', plistarray[1] class SkipThisVertex(Exception): pass def extractAntiSymmetricIndices(instring,funct) : terms = instring.strip(funct).strip(")").split(",") sign=1. for iy in range(0,len(terms)) : for ix in range(-1,-len(terms)+iy,-1) : swap = False if(len(terms[ix])==1 and len(terms[ix-1])==1) : swap = int(terms[ix])("\g<2>")', text) + text = open(os.path.join(model_dir, 'object_library.py'),'w').write(text) + text = open(os.path.join(model_dir, '__init__.py')).read() + mod = False + to_check = ['object_library', 'function_library'] + for lib in to_check: + if 'import %s' % lib in text: + continue + mod = True + text = "import %s \n" % lib + text + if mod: + open(os.path.join(model_dir, '__init__.py'),'w').write(text) + def convertFileToPython3(fName,names) : output="" inFile=open(fName) line=inFile.readline() isInit = "__init__.py" in fName tNames=[] for val in names : tNames.append(val) while line : # iteritems -> items line=line.replace("iteritems","items") - # fix imports + # fix imports if("import" in line) : for val in names : if("import %s" %val in line) : line=line.replace("import %s" %val, "from . import %s" %val) if(val in tNames) : tNames.remove(val) if("from %s" %val in line) : line=line.replace("from %s" %val, "from .%s" %val) # add brackets to print statements if("print" in line) : if line.strip()[0:5] == "print" : line=line.strip("\n").replace("print","print(")+")\n" output += line line=inFile.readline() inFile.close() if(isInit) : if "__init__" in tNames : tNames.remove("__init__") temp="" for val in tNames : temp+= "from . import %s\n" % val output = temp+output with open(fName,'w') as dest: dest.write(output)