Py8toJetInput Class Reference

#include <Py8toJetInput.h>

Inheritance diagram for Py8toJetInput:

Public Types
typedef Pythia8::Event	Event
typedef Pythia8::Particle	Particle

Public Member Functions
virtual const std::vector < fastjet::PseudoJet >	fillJetAlgoInput (const Event &, const Event &, const lhef::LHEEvent *lhee=0, const std::vector< int > *partonList=0)
	Py8toJetInput ()
void	setJetEtaMax (double max)
	~Py8toJetInput ()

Protected Types
enum	partonTypes { ID_TOP =6, ID_GLUON =21, ID_PHOTON =22 }

Protected Member Functions
int	getAncestor (int, const Event &, const Event &)

Protected Attributes
double	fJetEtaMax
std::vector< fastjet::PseudoJet >	fJetInput

Detailed Description

Definition at line 13 of file Py8toJetInput.h.

Member Typedef Documentation

typedef Pythia8::Event Py8toJetInput::Event

Definition at line 16 of file Py8toJetInput.h.

typedef Pythia8::Particle Py8toJetInput::Particle

Definition at line 17 of file Py8toJetInput.h.

Member Enumeration Documentation

enum Py8toJetInput::partonTypes

protected

Enumerator
ID_TOP
ID_GLUON
ID_PHOTON

Definition at line 29 of file Py8toJetInput.h.

{ ID_TOP=6, ID_GLUON=21, ID_PHOTON=22 };

Constructor & Destructor Documentation

Py8toJetInput::Py8toJetInput ( )

inline

Definition at line 19 of file Py8toJetInput.h.

: fJetEtaMax(10.) {}

Py8toJetInput::~Py8toJetInput ( )

inline

Definition at line 20 of file Py8toJetInput.h.

{}

Member Function Documentation

const std::vector< fastjet::PseudoJet > Py8toJetInput::fillJetAlgoInput ( const Event &  event, const Event &  workEvent, const lhef::LHEEvent *  lhee = 0, const std::vector< int > *  partonList = 0  )

virtual

Reimplemented in Py8toJetInputHEPEVT.

Definition at line 8 of file Py8toJetInput.cc.

References end, eta(), spr::find(), fJetEtaMax, fJetInput, lhef::LHEEvent::getHEPEUP(), i, ID_PHOTON, ID_TOP, j, lhef::HEPEUP::MOTHUP, pos, edm::shift, and ntuplemaker::status.

{

   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  )

protected

Definition at line 188 of file Py8toJetInput.cc.

References abs, pos, and ntuplemaker::status.

Referenced by Py8toJetInputHEPEVT::fillJetAlgoInput().

{

   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;

}

void Py8toJetInput::setJetEtaMax ( double  max )

inline

Definition at line 25 of file Py8toJetInput.h.

References fJetEtaMax, and max().

{ fJetEtaMax=max; return; }

Member Data Documentation

double Py8toJetInput::fJetEtaMax

protected

Definition at line 30 of file Py8toJetInput.h.

Referenced by fillJetAlgoInput(), and setJetEtaMax().

std::vector<fastjet::PseudoJet> Py8toJetInput::fJetInput

protected

Definition at line 34 of file Py8toJetInput.h.

Referenced by fillJetAlgoInput(), and Py8toJetInputHEPEVT::fillJetAlgoInput().