PowhegHooksBB4L.h

Go to the documentation of this file.

// PowhegHooksBB4L.h 
// Copyright (C) 2017 Silvia Ferrario Ravasio, Tomas Jezo, Paolo Nason, Markus Seidel
// inspired by PowhegHooks.h by Richard Corke
// adjusted to work with EmissionVetoHook1 in CMSSW by Alexander Grohsjean

#ifndef Pythia8_PowhegHooksBB4L_H
#define Pythia8_PowhegHooksBB4L_H

// Includes
#include "Pythia8/Pythia.h"
#include <cassert>
struct{
  int radtype;
} radtype_;

namespace Pythia8 {

  class PowhegHooksBB4L : public UserHooks {

  public:

    //--- Constructor and destructor -------------------------------------------
  PowhegHooksBB4L() : nFSRvetoBB4l(0) {}
    ~PowhegHooksBB4L() override {
      std::cout << "Number of FSR vetoed in BB4l = " << nFSRvetoBB4l << std::endl;
    }
    
    //--- Initialization -----------------------------------------------------------------------
    bool initAfterBeams() override {
      // settings of this class
      vetoFSREmission = settingsPtr->flag("POWHEG:bb4l:FSREmission:veto");
      onlyDistance1 = settingsPtr->flag("POWHEG:bb4l:FSREmission:onlyDistance1");
      dryRunFSR = settingsPtr->flag("POWHEG:bb4l:FSREmission:dryRun");
      vetoAtPL = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoAtPL");
      vetoQED = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoQED");
      vetoPartonLevel = settingsPtr->flag("POWHEG:bb4l:PartonLevel:veto");
      excludeFSRConflicting = settingsPtr->flag("POWHEG:bb4l:PartonLevel:excludeFSRConflicting");
      debug = settingsPtr->flag("POWHEG:bb4l:DEBUG");
      scaleResonanceVeto = settingsPtr->flag("POWHEG:bb4l:ScaleResonance:veto");
      vetoDipoleFrame = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoDipoleFrame");
      pTpythiaVeto = settingsPtr->flag("POWHEG:bb4l:FSREmission:pTpythiaVeto");
      //vetoProduction = (settingsPtr->mode("POWHEG:veto")==1);
      pTmin = settingsPtr->parm("POWHEG:bb4l:pTminVeto");
      return true;
    }


    //--- PROCESS LEVEL HOOK ---------------------------------------------------

    // called at the LHE level
    inline bool canVetoProcessLevel() override { return true; }
    inline bool doVetoProcessLevel(Event &e) override {
      
      // extract the radtype from the event comment
      stringstream ss;
      // use eventattribute as comments not filled when using edm input 
      //ss << infoPtr->getEventComments();
      ss << infoPtr->getEventAttribute("#rwgt");
      string temp;
      ss >> temp >> radtype_.radtype;
      assert (temp == "#rwgt");
      
      // find last top and the last anti-top in the record
      int i_top = -1, i_atop = -1;
      for (int i = 0; i < e.size(); i++) {
    if (e[i].id() == 6) i_top = i;
    if (e[i].id() == -6) i_atop = i;
      }
      if (i_top != -1) 
    topresscale = findresscale(i_top, e);
      else 
    topresscale = 1e30;
      if (i_top != -1) 
    atopresscale = findresscale(i_atop, e);
      else 
    atopresscale = 1e30;
      // initialize stuff
      doVetoFSRInit();
      // do not veto, ever
      return false;
    }

    //--- PARTON LEVEL HOOK ----------------------------------------------------

    // called after shower
    bool retryPartonLevel() override { return vetoPartonLevel || vetoAtPL; }
    inline bool canVetoPartonLevel() override { return vetoPartonLevel || vetoAtPL; }
    inline bool doVetoPartonLevel(const Event &e) override {
      if(radtype_.radtype==2)
    return false;
      if (debug){
    if (dryRunFSR && wouldVetoFsr) {
      double scale = getdechardness(vetoTopCharge, e);
      cout << "FSRdecScale = " << vetoDecScale << ", PLdecScale = " << scale << ", ratio = " << vetoDecScale/scale << endl;
    }
      }
      if (vetoPartonLevel) {
    double topdecscale = getdechardness(1, e);
    double atopdecscale = getdechardness(-1, e);
    if ((topdecscale > topresscale) || (atopdecscale > atopresscale)) {
      //if(dryRunFSR && ! wouldVetoFsr) mydatacontainer_.excludeEvent = excludeFSRConflicting?1:0;
      return true;
    }
    else 
      //if(dryRunFSR && wouldVetoFsr) mydatacontainer_.excludeEvent = excludeFSRConflicting?1:0;
      return false;
      }
      if (vetoAtPL) {
    if (dryRunFSR && wouldVetoFsr) return true;
    else return false;
      }
      return false;
    }

    //--- FSR EMISSION LEVEL HOOK ----------------------------------------------

    // FSR veto: this should be true if we want PowhegHooksBB4l veto in decay
    //           OR PowhegHooks veto in production. (The virtual method
    //           PowhegHooks::canVetoFSREmission has been replaced by
    //           PowhegHooksBB4L::canVetoFSREmission, so FSR veto in production
    //           must be handled here. ISR and MPI veto are instead still
    //           handled by PowhegHooks.)
    inline bool canVetoFSREmission() override { return vetoFSREmission; } // || vetoProduction; }
    inline bool doVetoFSREmission(int sizeOld, const Event &e, int iSys, bool inResonance) override {
      //VETO INSIDE THE RESONANCE //
      if (inResonance && vetoFSREmission) {
    int iRecAft = e.size() - 1;
    int iEmt = e.size() - 2;
    int iRadAft = e.size() - 3;
    int iRadBef = e[iEmt].mother1();

    // find the top resonance the radiator originates from
    int iTop = e[iRadBef].mother1();
    int distance = 1;
    while (abs(e[iTop].id()) != 6 && iTop > 0) {
      iTop = e[iTop].mother1();
      distance ++;
    }
    if (iTop == 0) {
      infoPtr->errorMsg("Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from top quark, not vetoing");
      return doVetoFSR(false,0,0);
      //return false;
    }
    int iTopCharge = (e[iTop].id()>0)?1:-1;

    // calculate the scale of the emission
    double scale;
    //using pythia pT definition ...
    if(pTpythiaVeto)
      scale = pTpythia(e, iRadAft, iEmt, iRecAft);
    //.. or using POWHEG pT definition
    else{
      Vec4 pr(e[iRadAft].p()), pe(e[iEmt].p()), pt(e[iTop].p()), prec(e[iRecAft].p()), psystem;
      // The computation of the POWHEG pT can be done in the top rest frame or in the diple one.
      // pdipole = pemt +prec +prad (after the emission)
      // For the first emission off the top resonance pdipole = pw +pb (before the emission) = ptop
      if(vetoDipoleFrame)
        psystem = pr+pe+prec;
      else
        psystem = pt;
            
      // gluon splitting into two partons
      if (e[iRadBef].id() == 21)
        scale = gSplittingScale(psystem, pr, pe);
      // quark emitting a gluon (or a photon)
      else if (abs(e[iRadBef].id()) == 5 && ((e[iEmt].id() == 21) && ! vetoQED) )
        scale = qSplittingScale(psystem, pr, pe);
      // other stuff (which we should not veto)
      else 
        scale = 0;
    }
        
    if (iTopCharge > 0) {
      if (onlyDistance1) {
        if ( debug && (distance == 1) && scale > topresscale && ! wouldVetoFsr)
          cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl;
        return doVetoFSR((distance == 1) && scale > topresscale,scale,iTopCharge);
      }
      else {
        if ( debug && scale > topresscale && ! wouldVetoFsr)
          cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl;
        return doVetoFSR(scale > topresscale,scale,iTopCharge);
      }
    }
    else if (iTopCharge < 0){
      if (onlyDistance1){
        if ( debug && (distance == 1) && scale > atopresscale && ! wouldVetoFsr)
          cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl;
        return doVetoFSR((distance == 1) && scale > atopresscale,scale,iTopCharge);
      }
      else {
        if ( debug && scale > topresscale && ! wouldVetoFsr)
          cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl;
        return doVetoFSR(scale > atopresscale,scale,iTopCharge);
      }
    }
    else {
      cout << "Bug in PohwgeHooksBB4l" << endl;     
    }
      }
      // VETO THE PRODUCTION PROCESS //
      // covered by multiuserhook, i.e. need to turn on EV1    
      // else if(!inResonance && vetoProduction){
      // return EmissionVetoHook1::doVetoFSREmission(sizeOld, e, iSys, inResonance);
      // }

      return false; 
    }

    inline bool doVetoFSR(bool condition, double scale, int iTopCharge)  {
      if(radtype_.radtype==2)
    return false;
      if (condition) {
    if (!wouldVetoFsr) {
      wouldVetoFsr = true;
      vetoDecScale = scale;
      vetoTopCharge = iTopCharge;
    }
    if (dryRunFSR) return false;
    else { 
      nFSRvetoBB4l++; 
      return true; 
    }
      }
      else return false;
    }

    inline void doVetoFSRInit() {
      wouldVetoFsr = false;
      vetoDecScale = -1;
      vetoTopCharge = 0;
    }

    //--- SCALE RESONANCE HOOK -------------------------------------------------
    // called before each resonance decay shower
    inline bool canSetResonanceScale() override { return scaleResonanceVeto; }
    // if the resonance is the (anti)top set the scale to:
    //  ---> (anti)top virtuality if radtype=2
    //  ---> (a)topresscale otherwise
    // if is not the top, set it to a big number
    inline double scaleResonance(int iRes, const Event &e) override {       
      if (e[iRes].id() == 6){
    if(radtype_.radtype == 2)
      return sqrt(e[iRes].m2Calc());
    else
      return topresscale;
      }
      else if (e[iRes].id() == -6){
    if(radtype_.radtype == 2)
      return sqrt(e[iRes].m2Calc());
    else
      return atopresscale;
      }
      else
    return pow(10.0,30.);
    }


    //--- Internal helper functions --------------------------------------------

    // Calculates the scale of the hardest emission from within the resonance system
    // translated by Markus Seidel modified by Tomas Jezo
    inline double findresscale( const int iRes, const Event& event) {
      double scale = 0.;

      int nDau = event[iRes].daughterList().size();

      if (nDau == 0) {
    // No resonance found, set scale to high value
    // Pythia will shower any MC generated resonance unrestricted
    scale = 1e30;
      }
      else if (nDau < 3) {
    // No radiating resonance found
    scale = pTmin;
      }
      else if (abs(event[iRes].id()) == 6) {
    // Find top daughters
    int idw = -1, idb = -1, idg = -1;
            
    for (int i = 0; i < nDau; i++) {
      int iDau = event[iRes].daughterList()[i];
      if (abs(event[iDau].id()) == 24) idw = iDau;
      if (abs(event[iDau].id()) ==  5) idb = iDau;
      if (abs(event[iDau].id()) == 21) idg = iDau;
    }
            
    // Get daughter 4-vectors in resonance frame
    Vec4 pw(event[idw].p());
    pw.bstback(event[iRes].p());
            
    Vec4 pb(event[idb].p());
    pb.bstback(event[iRes].p());
            
    Vec4 pg(event[idg].p());
    pg.bstback(event[iRes].p());
            
    // Calculate scale
    scale = sqrt(2*pg*pb*pg.e()/pb.e());
      }
      else {
    scale = 1e30;
      }

      return scale;
    }

    // The following routine will match daughters of particle `e[iparticle]` to an expected pattern specified via the list of expected particle PDG ID's `ids`, 
    // id wildcard is specified as 0 if match is obtained, the positions and the momenta of these particles are returned in vectors `positions` and `momenta` 
    // respectively
    // if exitOnExtraLegs==true, it will exit if the decay has more particles than expected, but not less
    inline bool match_decay(int iparticle, const Event &e, const vector<int> &ids, vector<int> &positions, vector<Vec4> &momenta, bool exitOnExtraLegs = true){
      // compare sizes
      if (e[iparticle].daughterList().size() != ids.size()) {
    if (exitOnExtraLegs && e[iparticle].daughterList().size() > ids.size()) {
      cout << "extra leg" << endl;
      exit(-1);
    }
    return false; 
      }
      // compare content
      for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) {
    int di = e[iparticle].daughterList()[i];
    if (ids[i] != 0 && e[di].id() != ids[i]) 
      return false;
      }
      // reset the positions and momenta vectors (because they may be reused)
      positions.clear();
      momenta.clear();
      // construct the array of momenta
      for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) {
    int di = e[iparticle].daughterList()[i];
    positions.push_back(di);
    momenta.push_back(e[di].p());
      }
      return true;
    }

    inline double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){
      p1.bstback(pt);
      p2.bstback(pt);
      return sqrt( 2*p1*p2*p2.e()/p1.e() );
    }

    inline double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){
      p1.bstback(pt);
      p2.bstback(pt);       
      return sqrt( 2*p1*p2*p1.e()*p2.e()/(pow(p1.e()+p2.e(),2)) );
    }

    // Routines to calculate the pT (according to pTdefMode) in a FS splitting:
    // i (radiator before) -> j (emitted after) k (radiator after)
    // For the Pythia pT definition, a recoiler (after) must be specified.
    // (INSPIRED BY pythia8F77_31.cc double pTpythia)
    inline double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch,
               int RecAfterBranch)
    {

      // Convenient shorthands for later
      Vec4 radVec = e[RadAfterBranch].p();
      Vec4 emtVec = e[EmtAfterBranch].p();
      Vec4 recVec = e[RecAfterBranch].p();
      int  radID  = e[RadAfterBranch].id();

      // Calculate virtuality of splitting
      Vec4 Q(radVec + emtVec);
      double Qsq = Q.m2Calc();
        

      // Mass term of radiator
      double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
    pow2(particleDataPtr->m0(radID)) : 0.;
    
      // z values for FSR 
      double z, pTnow;
      // Construct 2 -> 3 variables
      Vec4 sum = radVec + recVec + emtVec;
      double m2Dip = sum.m2Calc();
          
      double x1 = 2. * (sum * radVec) / m2Dip;
      double x3 = 2. * (sum * emtVec) / m2Dip;
      z     = x1 / (x1 + x3);
      pTnow = z * (1. - z);


      // Virtuality
      pTnow *= (Qsq - m2Rad);

      if (pTnow < 0.) {
    cout << "Warning: pTpythia was negative" << endl;
    return -1.;
      }
      else
    return(sqrt(pTnow));
    }

    inline double getdechardness(int topcharge, const Event &e){
      int tid = 6*topcharge, wid = 24*topcharge, bid = 5*topcharge, gid = 21, wildcard = 0;
      // find last top in the record
      int i_top = -1;
      Vec4 p_top, p_b, p_g, p_g1, p_g2;
      for (int i = 0; i < e.size(); i++) 
    if (e[i].id() == tid) {
      i_top = i;
      p_top = e[i].p();
    }
      if (i_top == -1) return -1.0;
                
      // summary of cases
      // 1.) t > W b
      //   a.) b > 3     ... error
      //   b.) b > b g   ... h = sqrt(2*p_g*p_b*p_g.e()/p_b.e())
      //   c.) b > other ... h = -1
      //   return h
      // 2.) t > W b g
      //   a.)   b > 3     ... error
      //   b.)   b > b g   ... h1 = sqrt(2*p_g*p_b*p_g.e()/p_b.e())
      //   c.)   b > other ... h1 = -1
      //   i.)   g > 3     ... error
      //   ii.)  g > 2     ... h2 = sqrt(2*p_g1*p_g2*p_g1.e()*p_g2.e()/(pow(p_g1.e(),2)+pow(p_g2.e(),2))) );
      //   iii.) g > other ... h2 = -1
      //   return max(h1,h2)
      // 3.) else ... error

      vector<Vec4> momenta;
      vector<int> positions;

      // 1.) t > b W
      if ( match_decay(i_top, e, vector<int> {wid, bid}, positions, momenta, false) ) {
    double h;
    int i_b = positions[1];
    // a.+b.) b > 3 or b > b g 
    if ( match_decay(i_b, e, vector<int> {bid, gid}, positions, momenta) )
      h = qSplittingScale(e[i_top].p(), momenta[0], momenta[1]);
    // c.) b > other
    else 
      h = -1;
    return h;
      } 
      // 2.) t > b W g
      else if ( match_decay(i_top, e, vector<int> {wid, bid, gid}, positions, momenta, false) ) {
    double h1, h2;
    int i_b = positions[1], i_g = positions[2];
    // a.+b.) b > 3 or b > b g
    if ( match_decay(i_b, e, vector<int> {bid, gid}, positions, momenta) )
      h1 = qSplittingScale(e[i_top].p(), momenta[0], momenta[1]);
    // c.) b > other
    else 
      h1 = -1;
    // i.+ii.) g > 3 or g > 2
    if ( match_decay(i_g, e, vector<int> {wildcard, wildcard}, positions, momenta) )
      h2 = gSplittingScale(e[i_top].p(), momenta[0], momenta[1]);
    // c.) b > other
    else 
      h2 = -1;
    return max(h1, h2);
      }
      // 3.) else
      else { 
    cout << "getdechardness" << endl;
    cout << "top at position " << i_top << endl;
    cout << "with " << e[i_top].daughterList().size() << " daughters " << endl;
    for (unsigned i = 0; i < e[i_top].daughterList().size(); i++) {
      int di = e[i_top].daughterList()[i];
      cout << "with daughter " << di << ": " << e[di].id() << endl;
    }
    exit(-1);
      }
    }



    //--------------------------------------------------------------------------

    // Functions to return information

    //  inline int    getNFSRveto() { return nFSRveto; } 

    //--------------------------------------------------------------------------

  private:
    // FSR emission veto flags
    bool vetoFSREmission, dryRunFSR, wouldVetoFsr, onlyDistance1, vetoAtPL, vetoQED;
    // Parton Level veto flags
    bool vetoPartonLevel, excludeFSRConflicting;
    // Scale Resonance veto flags
    double scaleResonanceVeto;
    // other flags
    bool debug;
    // internal: resonance scales
    double topresscale, atopresscale;
    // internal: inter veto communication
    double vetoDecScale;
    int vetoTopCharge;
    bool vetoDipoleFrame;
    bool pTpythiaVeto;
    //bool vetoProduction;
    double pTmin;
    // Statistics on vetos
    unsigned long int nFSRvetoBB4l; 

  };

  //==========================================================================

} // end namespace Pythia8

#endif // end Pythia8_PowhegHooksBB4L_H