PoWHEGVetoedShower.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
 
#include "UserHooksUtils.h"
#include "Pythia8_i/UserHooksFactory.h"
#include <stdexcept>
#include <iostream>
 
namespace Pythia8{
  class PoWHEGVetoedShower;
}
 
Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::PoWHEGVetoedShower> powhegVetoedShowerCreator("PoWHEGVetoedShower");
 
namespace Pythia8{
 
  class PoWHEGVetoedShower : public UserHooks{
 
  public:
 
    PoWHEGVetoedShower(): m_powhegScale(0.){
 
      std::cout<<"*******************************************************************"<<std::endl;
      std::cout<<"*                                                                 *"<<std::endl;
      std::cout<<"*  Using new PoWHEG author scale definition for QCD production!   *"<<std::endl;
      std::cout<<"*                                                                 *"<<std::endl;
      std::cout<<"*******************************************************************"<<std::endl;
      
    }
 
    ~PoWHEGVetoedShower(){}
 
    bool doVetoMPIStep(int nMPI, const Event &evt){
 
      // Only do anything on the first call.
      if(nMPI > 1) return false;
      
      m_powhegScale = infoPtr->QRen();;
      
      // momentum components to boost to CMS frame
      double pxCMS = 0.;
      double pyCMS = 0.;
      double pzCMS = 0.;
      double eCMS  = 0.;
      
      std::vector<Particle> powhegLegs;
      
      // Find the entries corresponding to outgoing legs from PoWHEG
      // start the loop at 1, since entry 0 represents the event as a whole
      for(int ii=1; ii != evt.size(); ++ii){
 
        // status -21 is the incoming hard partons
        if(evt[ii].status() == -21){
          pxCMS += evt[ii].px();
          pyCMS += evt[ii].py();
          pzCMS += evt[ii].pz();   
          eCMS  += evt[ii].e();
        }
        
        // here Event::isFinal refers to whether the particle can still decay/radiate, or whether it is an internal leg
        if(evt[ii].isFinal()){
          powhegLegs.push_back(Particle(evt[ii]));
        }
      }
      
      // compare the pT of each leg to the powheg scale.
      // Set the scale to the lowest (or leave the scale unchanged if it is already lower)
      for(std::vector<Particle>::const_iterator leg=powhegLegs.begin();
          leg != powhegLegs.end(); ++leg){
        if(leg->pT() < m_powhegScale )m_powhegScale = leg->pT();
      }
      
      double norm = -1./eCMS;
      pxCMS *= norm;
      pyCMS *= norm;
      pzCMS *= norm;
      
      // boost all outgoing legs to the CMS frame
      for(std::vector<Particle>::iterator leg=powhegLegs.begin();
          leg != powhegLegs.end(); ++leg){
        leg->bst(pxCMS, pyCMS, pzCMS);
      }
      
      for(std::vector<Particle>::const_iterator leg=powhegLegs.begin();
          leg != powhegLegs.end(); ++leg){
        
        if(leg->pT() < m_powhegScale )m_powhegScale = leg->pT();
        
        // calculate the pT relative to each other leg 
        // if any such pT is lower than the current scale, reset the scale to that value
        for(std::vector<Particle>::const_iterator otherLeg = powhegLegs.begin();
            otherLeg != powhegLegs.end(); ++otherLeg){
          if(otherLeg == leg) continue;
          
          double pTLeg = Pythia8_UserHooks::pTLeg(*leg, *otherLeg);
 
          if(pTLeg < m_powhegScale) m_powhegScale = pTLeg;
        }
      }
 
//      std::cout<<m_powhegScale<<"  "<<infoPtr->QRen()<<std::endl;
 
      
      return false;
    }
 
    bool doVetoISREmission(int, const Event &evt, int iSys){
 
      // only veto emissions from the hard system
      if(iSys != 0) return false;
      
      size_t emission = Pythia8_UserHooks::findLastISREmission(evt);
 
      // Veto emission above the veto scale
      if(evt[emission].pT() > m_powhegScale) return true;
      
      return false;
    }
    bool doVetoFSREmission(int, const Event &evt, int iSys, bool){
 
      // only veto emissions from the hard system
      if(iSys != 0) return false;
      
      size_t emission = Pythia8_UserHooks::findLastFSREmission(evt);
 
      // Veto if above the POWHEG scale
      if (evt[emission].pT() > m_powhegScale) return true;
      
      return false;
    }
 
    bool canVetoMPIStep(){return true;}
    int numberVetoMPIStep(){return 1;}
    bool canVetoISREmission(){return true;}
    bool canVetoFSREmission(){return true;}
 
  private:
 
    double m_powhegScale;
    
  };
}