LHAupAlpgen Class Reference

#include <GeneratorInput.h>

Inheritance diagram for LHAupAlpgen:

Public Member Functions
bool	fileFound ()
	LHAupAlpgen (const char *baseFNin, Pythia8::Info *infoPtrIn=NULL)
void	printParticles ()
bool	setEvent (int, double)
bool	setInit ()
	~LHAupAlpgen ()

Private Member Functions
bool	addResonances ()
bool	rescaleMomenta ()

Private Attributes
AlpgenPar	alpgenPar
string	baseFN
double	ebmupA
double	ebmupB
ifstream	ifsUnw
int	ihvy1
int	ihvy2
istream *	isUnw
int	lprup
double	mb
vector< Pythia8::LHAParticle >	myParticles
string	parFN
string	unwFN

Static Private Attributes
static const double	EWARNTHRESHOLD = 3e-3
static const double	INCOMINGMIN = 1e-3
static const bool	LHADEBUG = false
static const bool	LHADEBUGRESCALE = false
static const double	PTWARNTHRESHOLD = 1e-3
static const double	ZEROTHRESHOLD = 1e-10

Detailed Description

Definition at line 84 of file GeneratorInput.h.

Constructor & Destructor Documentation

LHAupAlpgen::LHAupAlpgen	(	const char *	baseFNin,
		Pythia8::Info *	infoPtrIn = NULL
	)

Definition at line 194 of file GeneratorInput.cc.

References alpgenPar, baseFN, ifsUnw, isUnw, NULL, parFN, AlpgenPar::parse(), edm::setPtr(), and unwFN.

  : baseFN(baseFNin), alpgenPar(infoPtrIn), isUnw(NULL) {

  // Set the info pointer if given
  if (infoPtrIn) setPtr(infoPtrIn);

  // Read in '_unw.par' file to get parameters
  ifstream  ifsPar;
  istream*  isPar = NULL;

  // Try gzip file first then normal file afterwards
#ifdef GZIPSUPPORT
  parFN = baseFN + "_unw.par.gz";
  isPar = openFile(parFN.c_str(), ifsPar);
  if (!ifsPar.is_open()) closeFile(isPar, ifsPar);
#endif
  if (isPar == NULL) {
    parFN = baseFN + "_unw.par";
    isPar = openFile(parFN.c_str(), ifsPar);
    if (!ifsPar.is_open()) {
      if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::LHAupAlpgen: "
          "cannot open parameter file", parFN);
      closeFile(isPar, ifsPar);
      return;
    }
  }

  // Read entire contents into string and close file
  string paramStr((istreambuf_iterator<char>(isPar->rdbuf())),
                   istreambuf_iterator<char>());

  // Make sure we reached EOF and not other error
  if (ifsPar.bad()) {
    if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::LHAupAlpgen: "
        "cannot read parameter file", parFN);
    return;
  }
  closeFile(isPar, ifsPar);

  // Parse file and set LHEF header
  alpgenPar.parse(paramStr);
  if (infoPtr) setInfoHeader("AlpgenPar", paramStr);

  // Open '.unw' events file (with possible gzip support)
#ifdef GZIPSUPPORT
  unwFN = baseFN + ".unw.gz";
  isUnw = openFile(unwFN.c_str(), ifsUnw);
  if (!ifsUnw.is_open()) closeFile(isUnw, ifsUnw);
#endif
  if (isUnw == NULL) {
    unwFN = baseFN + ".unw";
    isUnw = openFile(unwFN.c_str(), ifsUnw);
    if (!ifsUnw.is_open()) {
      if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::LHAupAlpgen: "
          "cannot open event file", unwFN);
      closeFile(isUnw, ifsUnw);
    }
  }
}

LHAupAlpgen::~LHAupAlpgen ( )

inline

Definition at line 90 of file GeneratorInput.h.

References ifsUnw, and isUnw.

{ closeFile(isUnw, ifsUnw); }

Member Function Documentation

bool LHAupAlpgen::addResonances ( )

private

Definition at line 484 of file GeneratorInput.cc.

References funct::abs(), i, customizeTrackingMonitorSeedNumber::idx, ihvy1, LHADEBUG, lprup, mb, myParticles, template_L1THistoryDQMService_cfg::outOfOrder, printParticles(), and mathSSE::sqrt().

Referenced by setEvent().

                                {

  // Temporary storage for resonance information
  int    idT, statusT, mother1T, mother2T, col1T, col2T;
  double pxT, pyT, pzT, eT, mT;
  // Leave tau and spin as default values
  double tauT = 0., spinT = 9.;

  // Alpgen process dependent parts. Processes:
  //    1 = W        + QQbar         + jets
  //    2 = Z/gamma* + QQbar         + jets
  //    3 = W                        + jets
  //    4 = Z/gamma*                 + jets
  //   10 = W        + c             + jets
  //   14 = W        + gamma         + jets
  //   15 = W        + QQbar + gamma + jets
  // When QQbar = ttbar, tops are not decayed in these processes.
  // Explicitly reconstruct W/Z resonances; assumption is that the
  // decay products are the last two particles.
  if (lprup <= 4 || lprup == 10 || lprup == 14 || lprup == 15) {
    // Decay products are the last two entries
    int i1 = myParticles.size() - 1, i2 = i1 - 1;

    // Follow 'alplib/alpsho.f' procedure to get ID
    if (myParticles[i1].idPart + myParticles[i2].idPart == 0)
      idT = 0;
    else
      idT = - (myParticles[i1].idPart % 2) - (myParticles[i2].idPart % 2);
    idT = (idT > 0) ? 24 : (idT < 0) ? -24 : 23;

    // Check that we get the expected resonance type; Z/gamma*
    if (lprup == 2 || lprup == 4) {
      if (idT != 23) {
        if (infoPtr) infoPtr->errorMsg("Error in "
            "LHAupAlpgen::addResonances: wrong resonance type in event");
        return false;
      }

    // W's
    } else {
      if (abs(idT) != 24) {
        if (infoPtr) infoPtr->errorMsg("Error in "
            "LHAupAlpgen::addResonances: wrong resonance type in event");
        return false;
      }
    }

    // Remaining input
    statusT  = 2;
    mother1T = 1;
    mother2T = 2;
    col1T = col2T = 0;
    pxT = myParticles[i1].pxPart + myParticles[i2].pxPart;
    pyT = myParticles[i1].pyPart + myParticles[i2].pyPart;
    pzT = myParticles[i1].pzPart + myParticles[i2].pzPart;
    eT  = myParticles[i1].ePart  + myParticles[i2].ePart;
    mT  = sqrt(eT*eT - pxT*pxT - pyT*pyT - pzT*pzT);
    myParticles.push_back(LHAParticle(
        idT, statusT, mother1T, mother2T, col1T, col2T,
        pxT, pyT, pzT, eT, mT, tauT, spinT, -1.));

    // Update decay product mothers (note array size as if from 1)
    myParticles[i1].mother1Part = myParticles[i2].mother1Part =
        myParticles.size();
    myParticles[i1].mother2Part = myParticles[i2].mother2Part = 0;

  // Processes:
  //    5 = nW + mZ + j gamma + lH + jets
  //    6 = QQbar         + jets    (QQbar = ttbar)
  //    8 = QQbar + Higgs + jets    (QQbar = ttbar)
  //   13 = top   + q               (topprc = 1)
  //   13 = top   + b               (topprc = 2)
  //   13 = top   + W     + jets    (topprc = 3)
  //   13 = top   + W     + b       (topprc = 4)
  //   16 = QQbar + gamma + jets    (QQbar = ttbar)
  //
  // When tops are present, they are decayed to Wb (both the W and b
  // are not given), with this W also decaying (decay products given).
  // The top is marked intermediate, the (intermediate) W is
  // reconstructed from its decay products, and the decay product mothers
  // updated. The final-state b is reconstructed from (top - W).
  //
  // W/Z resonances are given, as well as their decay products. The
  // W/Z is marked intermediate, and the decay product mothers updated.
  //
  // It is always assumed that decay products are at the end.
  // For processes 5 and 13, it is also assumed that the decay products
  // are in the same order as the resonances.
  // For processes 6, 8 and 16, the possibility of the decay products
  // being out-of-order is also taken into account.
  } else if ( ((lprup == 6 || lprup == 8 || lprup == 16) && ihvy1 == 6) ||
              lprup == 5 || lprup == 13) {

    // Go backwards through the particles looking for W/Z/top
    int idx = myParticles.size() - 1;
    for (int i = myParticles.size() - 1; i > -1; i--) {

      // W or Z
      if (myParticles[i].idPart == 23 ||
          abs(myParticles[i].idPart) == 24) {

        // Check that decay products and resonance match up
        int flav;
        if (myParticles[idx].idPart + myParticles[idx - 1].idPart == 0)
          flav = 0;
        else
          flav = - (myParticles[idx].idPart % 2)
                 - (myParticles[idx - 1].idPart % 2);
        flav = (flav > 0) ? 24 : (flav < 0) ? -24 : 23;
        if (flav != myParticles[i].idPart) {
          if (infoPtr)
            infoPtr->errorMsg("Error in LHAupAlpgen::addResonance: "
                "resonance does not match decay products");
          return false;
        }

        // Update status/mothers
        myParticles[i].statusPart      = 2;
        myParticles[idx  ].mother1Part = i + 1;
        myParticles[idx--].mother2Part = 0;
        myParticles[idx  ].mother1Part = i + 1;
        myParticles[idx--].mother2Part = 0;

      // Top
      } else if (abs(myParticles[i].idPart) == 6) {

        // Check that decay products and resonance match up
        int flav;
        if (myParticles[idx].idPart + myParticles[idx - 1].idPart == 0)
          flav = 0;
        else
          flav = - (myParticles[idx].idPart % 2)
                 - (myParticles[idx - 1].idPart % 2);
        flav = (flav > 0) ? 24 : (flav < 0) ? -24 : 23;

        bool outOfOrder = false, wrongFlavour = false;;
        if ( abs(flav) != 24 ||
             (flav ==  24 && myParticles[i].idPart !=  6) ||
             (flav == -24 && myParticles[i].idPart != -6) ) {

          // Processes 5 and 13, order should always be correct
          if (lprup == 5 || lprup == 13) {
            wrongFlavour = true;

          // Processes 6, 8 and 16, can have out of order decay products
          } else {

            // Go back two decay products and retry
            idx -= 2;
            if (myParticles[idx].idPart + myParticles[idx - 1].idPart == 0)
              flav = 0;
            else
              flav = - (myParticles[idx].idPart % 2)
                     - (myParticles[idx - 1].idPart % 2);
            flav = (flav > 0) ? 24 : (flav < 0) ? -24 : 23;

            // If still the wrong flavour then error
            if ( abs(flav) != 24 ||
                 (flav ==  24 && myParticles[i].idPart !=  6) ||
                 (flav == -24 && myParticles[i].idPart != -6) )
              wrongFlavour = true;
            else outOfOrder = true;
          }

          // Error if wrong flavour
          if (wrongFlavour) {
            if (infoPtr)
              infoPtr->errorMsg("Error in LHAupAlpgen::addResonance: "
                  "resonance does not match decay products");
            return false;
          }
        }

        // Mark t/tbar as now intermediate
        myParticles[i].statusPart = 2;

        // New intermediate W+/W-
        idT      = flav;
        statusT  = 2;
        mother1T = i + 1;
        mother2T = 0;
        col1T = col2T = 0;
        pxT = myParticles[idx].pxPart + myParticles[idx - 1].pxPart;
        pyT = myParticles[idx].pyPart + myParticles[idx - 1].pyPart;
        pzT = myParticles[idx].pzPart + myParticles[idx - 1].pzPart;
        eT  = myParticles[idx].ePart  + myParticles[idx - 1].ePart;
        mT  = sqrt(eT*eT - pxT*pxT - pyT*pyT - pzT*pzT);
        myParticles.push_back(LHAParticle(
            idT, statusT, mother1T, mother2T, col1T, col2T,
            pxT, pyT, pzT, eT, mT, tauT, spinT, -1.));

        // Update the decay product mothers
        myParticles[idx  ].mother1Part = myParticles.size();
        myParticles[idx--].mother2Part = 0;
        myParticles[idx  ].mother1Part = myParticles.size();
        myParticles[idx--].mother2Part = 0;

        // New final-state b/bbar
        idT     = (flav == 24) ? 5 : -5;
        statusT = 1;
        // Colour from top
        col1T   = myParticles[i].col1Part;
        col2T   = myParticles[i].col2Part;
        // Momentum from (t/tbar - W+/W-)
        pxT     = myParticles[i].pxPart - myParticles.back().pxPart;
        pyT     = myParticles[i].pyPart - myParticles.back().pyPart;
        pzT     = myParticles[i].pzPart - myParticles.back().pzPart;
        eT      = myParticles[i].ePart  - myParticles.back().ePart;
        mT      = sqrt(eT*eT - pxT*pxT - pyT*pyT - pzT*pzT);
        myParticles.push_back(LHAParticle(
            idT, statusT, mother1T, mother2T, col1T, col2T,
            pxT, pyT, pzT, eT, mT, tauT, spinT, -1.));

        // If decay products were out of order, reset idx to point
        // at correct decay products
        if (outOfOrder) idx += 4;

      } // if (abs(myParticles[i].idPart) == 6)
    } // for (i)


  // Processes:
  //    7 = QQbar + Q'Qbar' + jets (tops are not decayed)
  //    9 =                   jets
  //   11 = gamma           + jets
  //   12 = Higgs           + jets
  } else if (lprup == 7 || lprup == 9 || lprup == 11 || lprup == 12) {
    // Nothing to do for these processes
  }

  // For single top, set incoming b mass if necessary
  if (lprup == 13) for (int i = 0; i < 2; i++)
    if (abs(myParticles[i].idPart) == 5) {
      myParticles[i].mPart = mb;
      myParticles[i].ePart = sqrt(pow2(myParticles[i].pzPart) + pow2(mb));
    }

  // Debug output and done.
  if (LHADEBUG) printParticles();
  return true;

}

bool LHAupAlpgen::fileFound ( )

inline

Definition at line 93 of file GeneratorInput.h.

References isUnw, and NULL.

{ return (isUnw != NULL); }

void LHAupAlpgen::printParticles

(

)

Definition at line 457 of file GeneratorInput.cc.

References gather_cfg::cout, i, and myParticles.

Referenced by addResonances(), and rescaleMomenta().

                                 {

  cout << endl << "---- LHAupAlpgen particle listing begin ----" << endl;
  cout << scientific << setprecision(6);
  for (int i = 0; i < int(myParticles.size()); i++) {
    cout << setw(5)  << i
         << setw(5)  << myParticles[i].idPart
         << setw(5)  << myParticles[i].statusPart
         << setw(15) << myParticles[i].pxPart
         << setw(15) << myParticles[i].pyPart
         << setw(15) << myParticles[i].pzPart
         << setw(15) << myParticles[i].ePart
         << setw(15) << myParticles[i].mPart
         << setw(5)  << myParticles[i].mother1Part - 1
         << setw(5)  << myParticles[i].mother2Part - 1
         << setw(5)  << myParticles[i].col1Part
         << setw(5)  << myParticles[i].col2Part
         << endl;
  }
  cout << "----  LHAupAlpgen particle listing end  ----" << endl;
}

bool LHAupAlpgen::rescaleMomenta ( )

private

Definition at line 741 of file GeneratorInput.cc.

References a, funct::abs(), b, gather_cfg::cout, reco::dp, EWARNTHRESHOLD, i, j, LHADEBUGRESCALE, max(), myParticles, printParticles(), PTWARNTHRESHOLD, mathSSE::sqrt(), and ZEROTHRESHOLD.

Referenced by setEvent().

                                 {

  // Total momenta in/out
  int  nOut = 0;
  Vec4 pIn, pOut;
  for (int i = 0; i < int(myParticles.size()); i++) {
    Vec4 pNow = Vec4(myParticles[i].pxPart, myParticles[i].pyPart,
                     myParticles[i].pzPart, myParticles[i].ePart);
    if (i < 2) pIn += pNow;
    else if (myParticles[i].statusPart == 1) {
      nOut++;
      pOut += pNow;
    }
  }

  // pT out to match pT in. Split any imbalances over all outgoing
  // particles, and scale energies also to keep m^2 fixed.
  if (abs(pOut.pT() - pIn.pT()) > ZEROTHRESHOLD) {
    // Differences in px/py
    double pxDiff = (pOut.px() - pIn.px()) / nOut,
           pyDiff = (pOut.py() - pIn.py()) / nOut;

    // Warn if resulting changes above warning threshold
    if (pxDiff > PTWARNTHRESHOLD || pyDiff > PTWARNTHRESHOLD) {
      if (infoPtr) infoPtr->errorMsg("Warning in LHAupAlpgen::setEvent: "
          "large pT imbalance in incoming event");

      // Debug printout
      if (LHADEBUGRESCALE) {
        printParticles();
        cout << "pxDiff = " << pxDiff << ", pyDiff = " << pyDiff << endl;
      }
    }

    // Adjust all final-state outgoing
    pOut.reset();
    for (int i = 2; i < int(myParticles.size()); i++) {
      if (myParticles[i].statusPart != 1) continue;
      myParticles[i].pxPart -= pxDiff;
      myParticles[i].pyPart -= pyDiff;
      myParticles[i].ePart   = sqrt(max(0., pow2(myParticles[i].pxPart) +
          pow2(myParticles[i].pyPart) + pow2(myParticles[i].pzPart) +
          pow2(myParticles[i].mPart)));
      pOut += Vec4(myParticles[i].pxPart, myParticles[i].pyPart,
                   myParticles[i].pzPart, myParticles[i].ePart);
    }
  }

  // Differences in E/pZ and scaling factors
  double de = (pOut.e()  - pIn.e());
  double dp = (pOut.pz() - pIn.pz());
  double a  = 1 + (de + dp) / 2. / myParticles[0].ePart;
  double b  = 1 + (de - dp) / 2. / myParticles[1].ePart;

  // Warn if resulting energy changes above warning threshold.
  // Change in pz less than or equal to change in energy (incoming b
  // quark can have mass at this stage for process 13). Note that for
  // very small incoming momenta, the relative adjustment may be large,
  // but still small in absolute terms.
  if (abs(a - 1.) * myParticles[0].ePart > EWARNTHRESHOLD ||
      abs(b - 1.) * myParticles[1].ePart > EWARNTHRESHOLD) {
    if (infoPtr) infoPtr->errorMsg("Warning in LHAupAlpgen::setEvent: "
        "large rescaling factor");

    // Debug printout
    if (LHADEBUGRESCALE) {
      printParticles();
      cout << "de = " << de << ", dp = " << dp
           << ", a = " << a << ", b = " << b << endl
           << "Absolute energy change for incoming 0 = "
           << abs(a - 1.) * myParticles[0].ePart << endl
           << "Absolute energy change for incoming 1 = "
           << abs(b - 1.) * myParticles[1].ePart << endl;
    }
  }
  myParticles[0].ePart  *= a;
  myParticles[0].pzPart *= a;
  myParticles[1].ePart  *= b;
  myParticles[1].pzPart *= b;

  // Recalculate resonance four vectors
  for (int i = 0; i < int(myParticles.size()); i++) {
    if (myParticles[i].statusPart != 2) continue;

    // Only mothers stored in LHA, so go through all
    Vec4 resVec;
    for (int j = 0; j < int(myParticles.size()); j++) {
      if (myParticles[j].mother1Part - 1 != i) continue;
      resVec += Vec4(myParticles[j].pxPart, myParticles[j].pyPart,
                     myParticles[j].pzPart, myParticles[j].ePart);
    }

    myParticles[i].pxPart = resVec.px();
    myParticles[i].pyPart = resVec.py();
    myParticles[i].pzPart = resVec.pz();
    myParticles[i].ePart  = resVec.e();
  }

  return true;
}

bool LHAupAlpgen::setEvent	(	int	,
		double
	)

Definition at line 352 of file GeneratorInput.cc.

References funct::abs(), addResonances(), ebmupA, i, ifsUnw, INCOMINGMIN, isUnw, geometryCSVtoXML::line, lprup, max(), myParticles, nEvent, rescaleMomenta(), and mathSSE::sqrt().

                                      {

  // Read in the first line of the event
  int    nEvent, iProc, nParton;
  double Swgt, Sq;
  string line;
  if (!getline(*isUnw, line)) {
    // Read was bad
    if (ifsUnw.bad()) {
      if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::setEvent: "
          "could not read events from file");
      return false;
    }
    // End of file reached
    if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::setEvent: "
        "end of file reached");
    return false;
  }
  istringstream iss1(line);
  iss1 >> nEvent >> iProc >> nParton >> Swgt >> Sq;

  // Set the process details (ignore alphaQED and alphaQCD parameters)
  double wgtT = Swgt, scaleT = Sq;
  setProcess(lprup, wgtT, scaleT);

  // Incoming flavour and x information for later
  int    id1T, id2T;
  double x1T, x2T;
  // Temporary storage for read in parton information
  int    idT, statusT, mother1T, mother2T, col1T, col2T;
  double pxT, pyT, pzT, eT, mT;
  // Leave tau and spin as default values
  double tauT = 0., spinT = 9.;

  // Store particles locally at first so that resonances can be added
  myParticles.clear();

  // Now read in partons
  for (int i = 0; i < nParton; i++) {
    // Get the next line
    if (!getline(*isUnw, line)) {
      if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::setEvent: "
          "could not read events from file");
      return false;
    }
    istringstream iss2(line);

    // Incoming (flavour, colour, anticolour, pz)
    if (i < 2) {
      // Note that mothers will be set automatically by Pythia, and LHA
      // status -1 is for an incoming parton
      iss2 >> idT >> col1T >> col2T >> pzT;
      statusT  = -1;
      mother1T = mother2T = 0;
      pxT = pyT = mT = 0.;
      eT  = abs(pzT);

      // Adjust when zero pz/e
      if (pzT == 0.) {
        pzT = (i == 0) ? INCOMINGMIN : -INCOMINGMIN;
        eT  = INCOMINGMIN;
      }

    // Outgoing (flavour, colour, anticolour, px, py, pz, mass)
    } else {
      // Note that mothers 1 and 2 corresport to the incoming partons,
      // as set above, and LHA status +1 is for outgoing final state
      iss2 >> idT >> col1T >> col2T >> pxT >> pyT >> pzT >> mT;
      statusT  = 1;
      mother1T = 1;
      mother2T = 2;
      eT = sqrt(max(0., pxT*pxT + pyT*pyT + pzT*pzT + mT*mT));
    }

    // Add particle
    myParticles.push_back(LHAParticle(
        idT, statusT, mother1T, mother2T, col1T, col2T,
        pxT, pyT, pzT, eT, mT, tauT, spinT,-1.));
  }

  // Add resonances if required
  if (!addResonances()) return false;

  // Rescale momenta if required (must be done after full event
  // reconstruction in addResonances)
  if (!rescaleMomenta()) return false;

  // Pass particles on to Pythia
  for (size_t i = 0; i < myParticles.size(); i++)
    addParticle(myParticles[i]);

  // Set incoming flavour/x information and done
  id1T = myParticles[0].idPart;
  x1T  = myParticles[0].ePart / ebmupA;
  id2T = myParticles[1].idPart;
  x2T  = myParticles[1].ePart / ebmupA;
  setIdX(id1T, id2T, x1T, x2T);
  setPdf(id1T, id2T, x1T, x2T, 0., 0., 0., false);
  return true;
}

bool LHAupAlpgen::setInit

(

)

Definition at line 259 of file GeneratorInput.cc.

References alpgenPar, ebmupA, ebmupB, AlpgenPar::getParam(), AlpgenPar::getParamAsInt(), AlpgenPar::haveParam(), ihvy1, ihvy2, lprup, and mb.

                          {

  // Check that all required parameters are present
  if (!alpgenPar.haveParam("ih2") || !alpgenPar.haveParam("ebeam")  ||
      !alpgenPar.haveParam("hpc") || !alpgenPar.haveParam("xsecup") ||
      !alpgenPar.haveParam("xerrup")) {
    if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::setInit: "
        "missing input parameters");
    return false;
  }

  // Beam IDs
  int ih2 = alpgenPar.getParamAsInt("ih2");
  int idbmupA = 2212;
  int idbmupB = (ih2 == 1) ? 2212 : -2212;

  // Beam energies
  double ebeam = alpgenPar.getParam("ebeam");
  ebmupA = ebeam;
  ebmupB = ebmupA;

  // PDF group and set (at the moment not set)
  int pdfgupA = 0, pdfsupA = 0;
  int pdfgupB = 0, pdfsupB = 0;

  // Strategy is for unweighted events and xmaxup not important
  int    idwtup = 3;
  double xmaxup = 0.;

  // Get hard process code
  lprup = alpgenPar.getParamAsInt("hpc");

  // Check for unsupported processes
  if (lprup == 7 || lprup == 8 || lprup == 13) {
    if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::setInit: "
        "process not implemented");
    return false;
  }

  // Depending on the process code, get heavy flavour information:
  //    6 = QQbar           + jets
  //    7 = QQbar + Q'Qbar' + jets
  //    8 = QQbar + Higgs   + jets
  //   16 = QQbar + gamma   + jets
  if (lprup == 6 || lprup == 7 || lprup == 8 || lprup == 16) {
    if (!alpgenPar.haveParam("ihvy")) {
      if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::setInit: "
          "heavy flavour information not present");
      return false;
    }
    ihvy1 = alpgenPar.getParamAsInt("ihvy");

  } else ihvy1 = -1;
  if (lprup == 7) {
    if (!alpgenPar.haveParam("ihvy2")) {
      if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::setInit: "
          "heavy flavour information not present");
      return false;
    }
    ihvy2 = alpgenPar.getParamAsInt("ihvy2");
  } else ihvy2 = -1;
  // For single top (process 13), get b mass to set incoming
  mb = -1.;
  if (lprup == 13) {
    if (!alpgenPar.haveParam("mb")) {
      if (infoPtr) infoPtr->errorMsg("Error in LHAupAlpgen::setInit: "
          "heavy flavour information not present");
      return false;
    }
    mb = alpgenPar.getParam("mb");
  }

  // Set the beams
  setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
  setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
  setStrategy(idwtup);

  // Add the process
  double xsecup = alpgenPar.getParam("xsecup");
  double xerrup = alpgenPar.getParam("xerrup");
  addProcess(lprup, xsecup, xerrup, xmaxup);
  xSecSumSave = xsecup;
  xErrSumSave = xerrup;

  // All okay
  return true;
}

Member Data Documentation

AlpgenPar LHAupAlpgen::alpgenPar

private

Definition at line 112 of file GeneratorInput.h.

Referenced by LHAupAlpgen(), and setInit().

string LHAupAlpgen::baseFN

private

Definition at line 111 of file GeneratorInput.h.

Referenced by LHAupAlpgen().

double LHAupAlpgen::ebmupA

private

Definition at line 114 of file GeneratorInput.h.

Referenced by setEvent(), and setInit().

double LHAupAlpgen::ebmupB

private

Definition at line 114 of file GeneratorInput.h.

Referenced by setInit().

const double LHAupAlpgen::EWARNTHRESHOLD = 3e-3

staticprivate

Definition at line 123 of file GeneratorInput.h.

Referenced by rescaleMomenta().

ifstream LHAupAlpgen::ifsUnw

private

Definition at line 117 of file GeneratorInput.h.

Referenced by LHAupAlpgen(), setEvent(), and ~LHAupAlpgen().

int LHAupAlpgen::ihvy1

private

Definition at line 115 of file GeneratorInput.h.

Referenced by addResonances(), and setInit().

int LHAupAlpgen::ihvy2

private

Definition at line 115 of file GeneratorInput.h.

Referenced by setInit().

const double LHAupAlpgen::INCOMINGMIN = 1e-3

staticprivate

Definition at line 123 of file GeneratorInput.h.

Referenced by setEvent().

istream* LHAupAlpgen::isUnw

private

Definition at line 118 of file GeneratorInput.h.

Referenced by fileFound(), LHAupAlpgen(), setEvent(), and ~LHAupAlpgen().

const bool LHAupAlpgen::LHADEBUG = false

staticprivate

Definition at line 122 of file GeneratorInput.h.

Referenced by addResonances().

const bool LHAupAlpgen::LHADEBUGRESCALE = false

staticprivate

Definition at line 122 of file GeneratorInput.h.

Referenced by rescaleMomenta().

int LHAupAlpgen::lprup

private

Definition at line 113 of file GeneratorInput.h.

Referenced by addResonances(), setEvent(), and setInit().

double LHAupAlpgen::mb

private

Definition at line 116 of file GeneratorInput.h.

Referenced by addResonances(), and setInit().

vector<Pythia8::LHAParticle> LHAupAlpgen::myParticles

private

Definition at line 119 of file GeneratorInput.h.

Referenced by addResonances(), printParticles(), rescaleMomenta(), and setEvent().

string LHAupAlpgen::parFN

private

Definition at line 111 of file GeneratorInput.h.

Referenced by LHAupAlpgen().

const double LHAupAlpgen::PTWARNTHRESHOLD = 1e-3

staticprivate

Definition at line 123 of file GeneratorInput.h.

Referenced by rescaleMomenta().

string LHAupAlpgen::unwFN

private

Definition at line 111 of file GeneratorInput.h.

Referenced by LHAupAlpgen().

const double LHAupAlpgen::ZEROTHRESHOLD = 1e-10

staticprivate

Definition at line 123 of file GeneratorInput.h.

Referenced by rescaleMomenta().