#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.

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

Py8toJetInput::ID_GLUON

Definition: Py8toJetInput.h:29

Py8toJetInput::ID_PHOTON

Definition: Py8toJetInput.h:29

Py8toJetInput::ID_TOP

Definition: Py8toJetInput.h:29

Constructor & Destructor Documentation

Py8toJetInput::Py8toJetInput ( )

inline

Definition at line 19 of file Py8toJetInput.h.

19 : fJetEtaMax(10.) {}

Py8toJetInput::fJetEtaMax

double fJetEtaMax

Definition: Py8toJetInput.h:30

Py8toJetInput::~Py8toJetInput ( )

inline

Definition at line 20 of file Py8toJetInput.h.

20 {}

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, customizeTrackingMonitorSeedNumber::idx, j, lhef::HEPEUP::MOTHUP, edm::shift, and mps_update::status.

Referenced by JetMatchingHook::doVetoPartonLevel().

 {
 
    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 funct::abs(), counter, and mps_update::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;
 
 }