Py8toJetInput.cc

Go to the documentation of this file.

#include "GeneratorInterface/Pythia8Interface/plugins/Py8toJetInput.h"

#include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
#include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
#include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"

const std::vector<fastjet::PseudoJet> Py8toJetInput::fillJetAlgoInput(const Event& event,
                                                                      const Event& workEvent,
                                                                      const lhef::LHEEvent* lhee,
                                                                      const std::vector<int>* typeIdx) {
  fJetInput.clear();

  Event workEventJet = workEvent;

  const lhef::HEPEUP& hepeup = *lhee->getHEPEUP();

  std::set<int> typeSet[3];

  // FIXME !!!
  // This is not safe to assume it's 3 because we're passing in a pointer
  // and we do NOT know what the actuial size is. I'll have to improve it.
  //
  for (size_t i = 0; i < 3; i++) {
    typeSet[i].clear();
    for (size_t j = 0; j < typeIdx[i].size(); j++) {
      // HEPEUP and Py8 event record are shifted by 3
      // (system particle + 2 beam particles in Py8 event)
      // ... EXCEPT FOR THE DECAY PRODUCTS IF DONE AT THE ME LEVEL !!!!!
      // ... in such case, the decay productes get gutted out of the sequence
      // and get placed in Py8 Event later on...
      // so we need to figure out the shift
      int pos = typeIdx[i][j];
      int shift = 3;
      // skip the first 2 entrirs in HEPEUP - they're incoming partons
      for (int ip = 2; ip < pos; ip++) {
        // alternative can be: moth1 != 1 && moth2 !=2...
        // but moth1 == moth2 means pointer to the same mother t
        // that can only be a resonance, unless moth1==moth2==0
        //
        if (hepeup.MOTHUP[ip].first == hepeup.MOTHUP[ip].second) {
          shift--;
        }
      }
      pos += shift;
      // typeSet[i].insert( event[pos].daughter1() );
      typeSet[i].insert(pos);
    }
  }

  // int iTypeEnd = (typeIdx[1].empty()) ? 1 : 2;

  // --> FIXME !!!
  int iType = 0;     // only LIGHT jets for now
  int jetAllow = 0;  // hardcoded for now for the same value as is set in Steve's example
                     // at present, not even in use...
                     // int jetMatch = 0; // hardcoded for now for the same value as is set in Steve's example

  // Loop over particles and decide what to pass to the jet algorithm
  for (int i = 0; i < workEventJet.size(); ++i) {
    if (!workEventJet[i].isFinal())
      continue;

    // jetAllow option to disallow certain particle types
    if (jetAllow == 1) {
      // Original AG+Py6 algorithm explicitly excludes tops,
      // leptons and photons.
      int id = workEventJet[i].idAbs();
      if ((id >= 11 && id <= 16) || id == ID_TOP || id == ID_PHOTON) {
        workEventJet[i].statusNeg();
        continue;
      }
    }

    // Get the index of this particle in original event
    int idx = workEventJet[i].daughter1();

    // Start with particle idx, and afterwards track mothers
    while (true) {
      // Light jets
      if (iType == 0) {
        // Do not include if originates from heavy jet or 'other'
        if (typeSet[1].find(idx) != typeSet[1].end() || typeSet[2].find(idx) != typeSet[2].end()) {
          workEventJet[i].statusNeg();
          break;
        }

        // Made it to start of event record so done
        if (idx == 0) {
          break;
        }
        // Otherwise next mother and continue
        idx = event[idx].mother1();

        // Heavy jets
      } else if (iType == 1) {
        // Only include if originates from heavy jet
        if (typeSet[1].find(idx) != typeSet[1].end())
          break;

        // Made it to start of event record with no heavy jet mother,
        // so DO NOT include particle
        if (idx == 0) {
          workEventJet[i].statusNeg();
          break;
        }

        // Otherwise next mother and continue
        idx = event[idx].mother1();

      }  // if (iType)
    }    // while (true)
  }      // for (i)

  // For jetMatch = 2, insert ghost particles corresponding to
  // each hard parton in the original process
  /*
  if (jetMatch > 0) 
  {

    for (int i = 0; i < int(typeIdx[iType].size()); i++) 
    {
      // Get y/phi of the parton
      Vec4   pIn = eventProcess[typeIdx[iType][i]].p();
      double y   = Vec4y(pIn);
      double phi = pIn.phi();

      // Create a ghost particle and add to the workEventJet
      double e   = MG5_GHOSTENERGY;
      double e2y = exp(2. * y);
      double pz  = e * (e2y - 1.) / (e2y + 1.);
      double pt  = sqrt(e*e - pz*pz);
      double px  = pt * cos(phi);
      double py  = pt * sin(phi);
      workEventJet.append(Particle(ID_GLUON, 99, 0, 0, 0, 0, 0, 0,
                                px, py, pz, e, 0., 0, 9.));

    } // for (i)
  } // if (jetMatch == 2)
*/

  for (int i = 0; i < workEventJet.size(); i++) {
    // fisrt, weed out all entries marked with statusNeg();
    //
    if (workEventJet[i].status() < 0)
      continue;

    // now weed out all entries above etaMax
    // in principle, we can use etaJetMaxAlgo because it's set equal to etaJetMax
    // as for etaJetMax, it gets set to memain_.etaclmax
    //
    if (fabs(workEventJet[i].eta()) > fJetEtaMax)
      continue;

    // need to double check if native FastJet understands Py8 Event format
    // in general, PseudoJet gets formed from (px,py,pz,E)
    //
    fastjet::PseudoJet partTmp = workEventJet[i];
    fJetInput.push_back(partTmp);
  }

  return fJetInput;
}

int Py8toJetInput::getAncestor(int pos, const Event& fullEvent, const Event& workEvent) {
  int parentId = fullEvent[pos].mother1();
  int parentPrevId = 0;
  int counter = pos;

  while (parentId > 0) {
    if (parentId == fullEvent[counter].mother2())  // carbon copy, keep walking up
    {
      parentPrevId = parentId;
      counter = parentId;
      parentId = fullEvent[parentPrevId].mother1();
      continue;
    }

    // we get here if not a carbon copy

    // let's check if it's a normal process, etc.
    //
    if ((parentId < parentPrevId) || parentId < fullEvent[counter].mother2())  // normal process
    {
      // first of all, check if hard block
      if (abs(fullEvent[counter].status()) == 22 || abs(fullEvent[counter].status()) == 23) {
        // yes, it's the hard block
        // we got what we want, and can exit now !
        parentId = counter;
        break;
      } else {
        parentPrevId = parentId;
        parentId = fullEvent[parentPrevId].mother1();
      }
    } else if (parentId > parentPrevId ||
               parentId > pos)  // "circular"/"forward-pointing parent" - intermediate process
    {
      parentId = -1;
      break;
    }

    // additional checks... although we shouldn't be geeting here all that much...
    //
    if (abs(fullEvent[parentId].status()) == 22 || abs(fullEvent[parentId].status()) == 23)  // hard block
    {
      break;
    }
    if (abs(fullEvent[parentId].status()) < 22)  // incoming
    {
      parentId = -1;
      break;
    }
  }

  return parentId;
}

#include "HepMC/HEPEVT_Wrapper.h"
#include <cassert>

const std::vector<fastjet::PseudoJet> Py8toJetInputHEPEVT::fillJetAlgoInput(const Event& event,
                                                                            const Event& workEvent,
                                                                            const lhef::LHEEvent* lhee,
                                                                            const std::vector<int>*) {
  fJetInput.clear();

  HepMC::HEPEVT_Wrapper::zero_everything();

  // service container for further mother-daughters links
  //
  std::vector<int> Py8PartonIdx;  // position of original (LHE) partons in Py8::Event
  Py8PartonIdx.clear();
  std::vector<int> HEPEVTPartonIdx;  // position of LHE partons in HEPEVT (incl. ME-generated decays)
  HEPEVTPartonIdx.clear();

  // general counter
  //
  int index = 0;

  int Py8PartonCounter = 0;
  int HEPEVTPartonCounter = 0;

  // find the fisrt parton that comes from LHE (ME-generated)
  // skip the incoming particles/partons
  for (int iprt = 1; iprt < event.size(); iprt++) {
    const Particle& part = event[iprt];
    if (abs(part.status()) < 22)
      continue;  // below 10 is "service"
                 // 11-19 are beam particles; below 10 is "service"
                 // 21 is incoming partons
    Py8PartonCounter = iprt;
    break;
  }

  const lhef::HEPEUP& hepeup = *lhee->getHEPEUP();
  // start the counter from 2, because we don't want the incoming particles/oartons !
  for (int iprt = 2; iprt < hepeup.NUP; iprt++) {
    index++;
    HepMC::HEPEVT_Wrapper::set_id(index, hepeup.IDUP[iprt]);
    HepMC::HEPEVT_Wrapper::set_status(index, 2);
    HepMC::HEPEVT_Wrapper::set_momentum(
        index, hepeup.PUP[iprt][0], hepeup.PUP[iprt][1], hepeup.PUP[iprt][2], hepeup.PUP[iprt][4]);
    HepMC::HEPEVT_Wrapper::set_mass(index, hepeup.PUP[iprt][4]);
    // --> FIXME HepMC::HEPEVT_Wrapper::set_position( index, part.xProd(), part.yProd(), part.zProd(), part.tProd() );
    HepMC::HEPEVT_Wrapper::set_parents(index, 0, 0);  // NO, not anymore to the "system particle"
    HepMC::HEPEVT_Wrapper::set_children(index, 0, 0);
    if (hepeup.MOTHUP[iprt].first > 2 &&
        hepeup.MOTHUP[iprt].second > 2)  // decay from LHE, will NOT show at the start of Py8 event !!!
    {
      HEPEVTPartonCounter++;
      continue;
    }
    Py8PartonIdx.push_back(Py8PartonCounter);
    Py8PartonCounter++;
    HEPEVTPartonIdx.push_back(HEPEVTPartonCounter);
    HEPEVTPartonCounter++;
  }

  HepMC::HEPEVT_Wrapper::set_number_entries(index);

  // now that the initial partons are in, attach parton-level from Pythia8
  // do NOT reset index as we need to *add* more particles sequentially
  //
  for (int iprt = 1; iprt < workEvent.size(); iprt++)  // from 0-entry (system) or from 1st entry ???
  {
    const Particle& part = workEvent[iprt];

    //      if ( part.status() != 62 ) continue;
    if (part.status() < 51)
      continue;
    index++;
    HepMC::HEPEVT_Wrapper::set_id(index, part.id());

    // HepMC::HEPEVT_Wrapper::set_status( index, event.statusHepMC(iprt) );
    HepMC::HEPEVT_Wrapper::set_status(index, 1);
    HepMC::HEPEVT_Wrapper::set_momentum(index, part.px(), part.py(), part.pz(), part.e());
    HepMC::HEPEVT_Wrapper::set_mass(index, part.m());
    HepMC::HEPEVT_Wrapper::set_position(index, part.xProd(), part.yProd(), part.zProd(), part.tProd());
    HepMC::HEPEVT_Wrapper::set_parents(index, 0, 0);  // just set to 0 like in Py6...
                                                      // although for some, mother will need to be re-set properly !
    HepMC::HEPEVT_Wrapper::set_children(index, 0, 0);

    // now refine mother-daughters links, where applicable

    int parentId = getAncestor(part.daughter1(), event, workEvent);

    if (parentId <= 0)
      continue;

    for (int idx = 0; idx < (int)Py8PartonIdx.size(); idx++) {
      if (parentId == Py8PartonIdx[idx]) {
        int idx1 = HEPEVTPartonIdx[idx];
        HepMC::HEPEVT_Wrapper::set_parents(index, idx1 + 1, idx1 + 1);
        break;
      }
    }
  }

  HepMC::HEPEVT_Wrapper::set_number_entries(index);

  // --> FIXME   HepMC::HEPEVT_Wrapper::set_event_number( fEventNumber ); // well, if you know it... well, it's one of the counters...

  //   HepMC::HEPEVT_Wrapper::print_hepevt();

  return fJetInput;
}