Compatibility with pythia8 212 and features from 76X

FastSimulation/ParticleDecay/plugins/TestPythiaDecays.cc

@@ -76,7 +76,7 @@ class TestPythiaDecays : public edm::stream::EDAnalyzer <> {
   std::map<int,TH1D*> h_br;
   std::map<int,TH1D*> h_br_ref;
 
-  std::map<int,std::vector<string> > knownDecayModes;
+  std::map<int,std::vector<std::string> > knownDecayModes;
 
   Pythia8::Pythia * pythia;
   std::string outputFile;
@@ -182,7 +182,7 @@ TestPythiaDecays::TestPythiaDecays(const edm::ParameterSet& iConfig)
     h_br[pid]->SetCanExtend(TH1::kAllAxes);
     h_br_ref[pid] = (TH1D*)(h_br[pid]->Clone(strstr.str().c_str()));
     h_br_ref[pid]->SetTitle(h_br_ref[pid]->GetName());
-    knownDecayModes[pid] = vector<string>();
+    knownDecayModes[pid] = std::vector<std::string>();
     for(int d = 0;d<pd->sizeChannels();++d){
       Pythia8::DecayChannel & channel = pd->channel(d);
       std::vector<int> prod;
@@ -334,7 +334,7 @@ TestPythiaDecays::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetu
      size_t parentIndex = it->first;
      const SimTrack & parent = simtracks->at(parentIndex);
      int pid = abs(parent.type());
-     vector<size_t> & childIndices = it->second;
+     std::vector<size_t> & childIndices = it->second;
      if(childIndices.size() == 0)
        continue;
 

GeneratorInterface/Pythia8Interface/BuildFile.xml

@@ -6,6 +6,7 @@
 <use name="GeneratorInterface/Core"/>
 <use name="boost"/>
 <use name="pythia8"/>
+<use name="evtgen"/>
 <use name="hepmc"/>
 <use name="clhep"/>
 <export>

GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h

@@ -12,6 +12,8 @@
 #include <Pythia8/Pythia.h>
 #include <Pythia8Plugins/HepMC2.h>
 
+class EvtGenDecays;
+
 namespace CLHEP {
   class HepRandomEngine;
 }
@@ -51,6 +53,13 @@ namespace gen {
 	 unsigned int                   maxEventsToPrint;
          HepMC::IO_AsciiParticles*      ascii_io;
 
+         // EvtGen plugin
+         //
+         bool useEvtGen;
+         EvtGenDecays* evtgenDecays;
+         std::string evtgenDecFile;
+         std::string evtgenPdlFile;
+
       private:
 
          P8RndmEngine p8RndmEngine_;

GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.cc

@@ -1,9 +1,13 @@
+#include "Pythia8/Pythia.h"
+
+using namespace Pythia8;
+
 #include "GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h"
 
 #include "FWCore/ServiceRegistry/interface/Service.h"
-void EmissionVetoHook1::fatalEmissionVeto(string message) {
+void EmissionVetoHook1::fatalEmissionVeto(std::string message) {
   throw edm::Exception(edm::errors::Configuration,"Pythia8Interface")
-    << "EmissionVeto: " << message << endl;
+    << "EmissionVeto: " << message << std::endl;
 }
 
 //--------------------------------------------------------------------------
@@ -56,14 +60,14 @@ double EmissionVetoHook1::pTpythia(const Pythia8::Event &e, int RadAfterBranch,
 
   // Can get negative pT for massive splittings
   if (pTnow < 0.) {
-    cout << "Warning: pTpythia was negative" << endl;
+    std::cout << "Warning: pTpythia was negative" << std::endl;
     return -1.;
   }
 
 #ifdef DBGOUTPUT
-  cout << "pTpythia: rad = " << RadAfterBranch << ", emt = "
+  std::cout << "pTpythia: rad = " << RadAfterBranch << ", emt = "
        << EmtAfterBranch << ", rec = " << RecAfterBranch
-       << ", pTnow = " << sqrt(pTnow) << endl;
+       << ", pTnow = " << sqrt(pTnow) << std::endl;
 #endif
 
   // Return pT
@@ -97,13 +101,13 @@ double EmissionVetoHook1::pTpowheg(const Pythia8::Event &e, int i, int j, bool F
 
   // Check result
   if (pTnow < 0.) {
-    cout << "Warning: pTpowheg was negative" << endl;
+    std::cout << "Warning: pTpowheg was negative" << std::endl;
     return -1.;
   }
 
 #ifdef DBGOUTPUT
-  cout << "pTpowheg: i = " << i << ", j = " << j
-       << ", pTnow = " << pTnow << endl;
+  std::cout << "pTpowheg: i = " << i << ", j = " << j
+       << ", pTnow = " << pTnow << std::endl;
 #endif
 
    return pTnow;
@@ -163,7 +167,7 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i
           // ISR - only done once as just kinematical pT
           if (!FSR) {
             pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
-            if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
+            if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
 
           // FSR - try all outgoing partons from system before branching 
           // as i. Note that for the hard system, there is no 
@@ -182,7 +186,7 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i
               pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0) 
                 ? false : FSR);
               if (pTnow > 0.) pTemt = (pTemt < 0) 
-                ? pTnow : min(pTemt, pTnow);
+                ? pTnow : std::min(pTemt, pTnow);
             } // for (iMem)
 
           } // if (!FSR)
@@ -193,9 +197,9 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i
           // ISR - other incoming as recoiler
           if (!FSR) {
             pTnow = pTpythia(e, iInA, jNow, iInB, FSR);
-            if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
+            if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
             pTnow = pTpythia(e, iInB, jNow, iInA, FSR);
-            if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
+            if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
 
           // FSR - try all final-state coloured partons as radiator
           //       after emission (k).
@@ -208,18 +212,18 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i
               // Try two incoming.
               pTnow = pTpythia(e, kNow, jNow, iInA, FSR);
               if (pTnow > 0.) pTemt = (pTemt < 0) 
-                ? pTnow : min(pTemt, pTnow);
+                ? pTnow : std::min(pTemt, pTnow);
               pTnow = pTpythia(e, kNow, jNow, iInB, FSR);
               if (pTnow > 0.) pTemt = (pTemt < 0) 
-                ? pTnow : min(pTemt, pTnow);
+                ? pTnow : std::min(pTemt, pTnow);
 
               // Try all other outgoing.
               for (int rNow = 0; rNow < e.size(); rNow++) {
                 if (rNow == kNow || rNow == jNow ||
                     !e[rNow].isFinal() || e[rNow].colType() == 0) continue;
                 pTnow = pTpythia(e, kNow, jNow, rNow, FSR);
                 if (pTnow > 0.) pTemt = (pTemt < 0) 
-                  ? pTnow : min(pTemt, pTnow);
+                  ? pTnow : std::min(pTemt, pTnow);
               } // for (rNow)
 
             } // for (kNow)
@@ -230,9 +234,9 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i
   } // for (xSR)
 
 #ifdef DBGOUTPUT
-  cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k
+  std::cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k
        << ", r = " << r << ", xSR = " << xSRin
-       << ", pTemt = " << pTemt << endl;
+       << ", pTemt = " << pTemt << std::endl;
 #endif
 
   return pTemt;
@@ -274,7 +278,7 @@ bool EmissionVetoHook1::doVetoMPIStep(int nMPI, const Pythia8::Event &e) {
       first = ip;
       break;
     }
-    if(first < 0) fatalEmissionVeto(string("signal particles not found"));
+    if(first < 0) fatalEmissionVeto(std::string("signal particles not found"));
     for(int ip = first; ip < e.size(); ip++) {
       myid = e[ip].id();
       if(abs(myid) < 6 || abs(myid) == 21) continue;
@@ -285,7 +289,7 @@ bool EmissionVetoHook1::doVetoMPIStep(int nMPI, const Pythia8::Event &e) {
 
   // Extra check that we have the correct final state
   if (count != nFinal && count != nFinal + 1)
-    fatalEmissionVeto(string("Wrong number of final state particles in event"));
+    fatalEmissionVeto(std::string("Wrong number of final state particles in event"));
 
   // Flag if POWHEG radiation present and index
   bool isEmt = (count == nFinal) ? false : true;
@@ -314,9 +318,9 @@ bool EmissionVetoHook1::doVetoMPIStep(int nMPI, const Pythia8::Event &e) {
   }
 
   if(Verbosity)
-    cout << "doVetoMPIStep: QFac = " << infoPtr->QFac()
+    std::cout << "doVetoMPIStep: QFac = " << infoPtr->QFac()
          << ", QRen = " << infoPtr->QRen()
-         << ", pThard = " << pThard << endl << endl;
+         << ", pThard = " << pThard << std::endl << std::endl;
 
   // Initialise other variables
   accepted   = false;
@@ -347,12 +351,12 @@ bool EmissionVetoHook1::doVetoISREmission(int, const Pythia8::Event &e, int iSys
   }
   if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
     e.list();
-    fatalEmissionVeto(string("Couldn't find Pythia ISR emission"));
+    fatalEmissionVeto(std::string("Couldn't find Pythia ISR emission"));
   }
 
   // pTemtMode == 0: pT of emitted w.r.t. radiator
-  // pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing)
-  // pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing)
+  // pTemtMode == 1: std::min(pT of emitted w.r.t. all incoming/outgoing)
+  // pTemtMode == 2: std::min(pT of all outgoing w.r.t. all incoming/outgoing)
   int xSR      = (pTemtMode == 0) ? 0       : -1;
   int i        = (pTemtMode == 0) ? iRadAft : -1;
   int j        = (pTemtMode != 2) ? iEmt    : -1;
@@ -361,7 +365,7 @@ bool EmissionVetoHook1::doVetoISREmission(int, const Pythia8::Event &e, int iSys
   double pTemt = pTcalc(e, i, j, k, r, xSR);
 
 #ifdef DBGOUTPUT
-  cout << "doVetoISREmission: pTemt = " << pTemt << endl << endl;
+  std::cout << "doVetoISREmission: pTemt = " << pTemt << std::endl << std::endl;
 #endif
 
   // Veto if pTemt > pThard
@@ -396,13 +400,13 @@ bool EmissionVetoHook1::doVetoFSREmission(int, const Pythia8::Event &e, int iSys
   if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) ||
        e[iEmt].status() != 51 || e[iRadAft].status() != 51) {
     e.list();
-    fatalEmissionVeto(string("Couldn't find Pythia FSR emission"));
+    fatalEmissionVeto(std::string("Couldn't find Pythia FSR emission"));
   }
 
   // Behaviour based on pTemtMode:
   //  0 - pT of emitted w.r.t. radiator before
-  //  1 - min(pT of emitted w.r.t. all incoming/outgoing)
-  //  2 - min(pT of all outgoing w.r.t. all incoming/outgoing)
+  //  1 - std::min(pT of emitted w.r.t. all incoming/outgoing)
+  //  2 - std::min(pT of all outgoing w.r.t. all incoming/outgoing)
   int xSR = (pTemtMode == 0) ? 1       : -1;
   int i   = (pTemtMode == 0) ? iRadBef : -1;
   int k   = (pTemtMode == 0) ? iRadAft : -1;
@@ -421,11 +425,11 @@ bool EmissionVetoHook1::doVetoFSREmission(int, const Pythia8::Event &e, int iSys
 
     for (int jLoop = 0; jLoop < 2; jLoop++) {
       if      (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR);
-      else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR));
+      else if (jLoop == 1) pTemt = std::min(pTemt, pTcalc(e, i, j, k, r, xSR));
 
       // For emittedMode == 3, have tried iRadAft, now try iEmt
       if (emittedMode != 3) break;
-      if (k != -1) swap(j, k); else j = iEmt;
+      if (k != -1) std::swap(j, k); else j = iEmt;
     }
 
   // If pTemtMode is 2, then try all final-state partons as emitted
@@ -435,7 +439,7 @@ bool EmissionVetoHook1::doVetoFSREmission(int, const Pythia8::Event &e, int iSys
   }
 
 #ifdef DBGOUTPUT
-  cout << "doVetoFSREmission: pTemt = " << pTemt << endl << endl;
+  std::cout << "doVetoFSREmission: pTemt = " << pTemt << std::endl << std::endl;
 #endif
 
   // Veto if pTemt > pThard
@@ -459,8 +463,8 @@ bool EmissionVetoHook1::doVetoMPIEmission(int, const Pythia8::Event &e) {
   if (MPIvetoOn) {
     if (e[e.size() - 1].pT() > pTMPI) {
 #ifdef DBGOUTPUT
-      cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT()
-           << ", pTMPI = " << pTMPI << endl << endl;
+      std::cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT()
+           << ", pTMPI = " << pTMPI << std::endl << std::endl;
 #endif
       return true;
     }

GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h

@@ -1,5 +1,3 @@
-#include "Pythia8/Pythia.h"
-
 class EmissionVetoHook1 : public Pythia8::UserHooks {
 
 public:  
@@ -15,8 +13,8 @@ class EmissionVetoHook1 : public Pythia8::UserHooks {
                     MPIvetoOn(MPIvetoOnIn), nISRveto(0), nFSRveto(0),
                     Verbosity(VerbosityIn) {}
  ~EmissionVetoHook1() {
-    cout << "Number of ISR vetoed = " << nISRveto << endl;
-    cout << "Number of FSR vetoed = " << nFSRveto << endl;
+    std::cout << "Number of ISR vetoed = " << nISRveto << std::endl;
+    std::cout << "Number of FSR vetoed = " << nFSRveto << std::endl;
   }
 
 //--------------------------------------------------------------------------
@@ -34,7 +32,7 @@ class EmissionVetoHook1 : public Pythia8::UserHooks {
   bool canVetoMPIEmission() { return MPIvetoOn; }
   bool doVetoMPIEmission(int, const Pythia8::Event &e);
 
-  void fatalEmissionVeto(string message);  
+  void fatalEmissionVeto(std::string message);
 
   double pTpythia(const Pythia8::Event &e, int RadAfterBranch,
                   int EmtAfterBranch, int RecAfterBranch, bool FSR);

GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.cc

@@ -42,10 +42,13 @@ bool LHAupLesHouches::setInit()
     infoPtr->setHeader("slha",slhaheader);
   }  
 
-  //work around missing initialization inside pythia8
-  infoPtr->eventAttributes = new std::map<std::string, std::string >;
-
-
+  //will be used to work around missing initialization inside pythia8
+  if(!fEvAttributes) {
+    fEvAttributes = new std::map<std::string, std::string >;
+  } else {
+    fEvAttributes->clear();
+  }
+
   return true;
 }
 
@@ -91,8 +94,13 @@ bool LHAupLesHouches::setEvent(int inProcId)
                 hepeup.SPINUP[i],scalein);
   }
 
-  infoPtr->eventAttributes->clear();
-
+  if(!infoPtr->eventAttributes) {
+    fEvAttributes->clear();
+    infoPtr->eventAttributes = fEvAttributes;
+  } else {
+    infoPtr->eventAttributes->clear();
+  }
+
   //fill parton multiplicities if available
   int npLO = event->npLO();
   int npNLO = event->npNLO();

GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.h

@@ -21,7 +21,7 @@
 
 class LHAupLesHouches : public Pythia8::LHAup {
   public:
-    LHAupLesHouches() : setScalesFromLHEF_(false) {;}
+    LHAupLesHouches() : setScalesFromLHEF_(false),fEvAttributes(0) {;}
 
     //void loadRunInfo(const boost::shared_ptr<lhef::LHERunInfo> &runInfo)
     void loadRunInfo(lhef::LHERunInfo* runInfo)
@@ -33,6 +33,8 @@ class LHAupLesHouches : public Pythia8::LHAup {
 
     void setScalesFromLHEF(bool b) { setScalesFromLHEF_ = b; }
 
+    ~LHAupLesHouches() {if(fEvAttributes) delete fEvAttributes;}
+
   private:
 
     bool setInit();
@@ -46,4 +48,5 @@ class LHAupLesHouches : public Pythia8::LHAup {
     // Flag to set particle production scales or not.
     bool setScalesFromLHEF_;
 
+    std::map<std::string, std::string> * fEvAttributes;
 };

GeneratorInterface/Pythia8Interface/plugins/Py8EGun.cc

@@ -4,6 +4,10 @@
 
 #include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h"
 
+// EvtGen plugin
+//
+#include "Pythia8Plugins/EvtGen.h"
+
 namespace gen {
 
 class Py8EGun : public Py8GunBase {
@@ -80,18 +84,20 @@ bool Py8EGun::generatePartonsAndHadronize()
       {
          if ( (fMasterGen->particleData).isParticle( -particleID ) )
 	 {
-	    (fMasterGen->event).append( -particleID, 1, 0, 0, px, py, pz, ee, mass );
+	    (fMasterGen->event).append( -particleID, 1, 0, 0, -px, -py, -pz, ee, mass );
 	 }
 	 else
 	 {
-	    (fMasterGen->event).append( particleID, 1, 0, 0, px, py, pz, ee, mass );
+	    (fMasterGen->event).append( particleID, 1, 0, 0, -px, -py, -pz, ee, mass );
 	 }
       }
 
    }
 
    if ( !fMasterGen->next() ) return false;
 
+   if (evtgenDecays) evtgenDecays->decay();
+
    event().reset(new HepMC::GenEvent);
    return toHepMC.fill_next_event( fMasterGen->event, event().get() );