WprimeFlat.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
 
#include "UserSetting.h"
#include "Pythia8_i/UserHooksFactory.h"
#include "Pythia8/PhaseSpace.h"
 
#include <stdexcept>
 
namespace Pythia8{
  class WprimeFlat;
}
 
Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::WprimeFlat> WprimeFlatCreator("WprimeFlat");
 
namespace Pythia8 {
 
  class WprimeFlat : public UserHooks {
    
  public:
    
    // Constructor.
    WprimeFlat(): m_energyMode("WprimeFlat:EnergyMode", 8){}
    
    // Destructor.
    ~WprimeFlat(){}
    
    // Allow process cross section to be modified...
    virtual bool canModifySigma() {return true;}
  
    // ...which gives access to the event at the trial level, before selection.
    virtual double multiplySigmaBy(const SigmaProcess* sigmaProcessPtr, 
                   const PhaseSpace* phaseSpacePtr, 
                   bool /* inEvent */) {
      // All events should be 2 -> 1, but kill them if not.
      if (sigmaProcessPtr->nFinal() != 1) return 0.; 
      
      // Weight cross section with BW propagator, i.e. to remove it.
      // (inEvent = false for initialization).
      // No inEvent criteria, want weight both for cross section 
      // and MC generation.
      int idRes   = sigmaProcessPtr->resonanceA();
      double mRes = particleDataPtr->m0(idRes);
      double wRes = particleDataPtr->mWidth(idRes);
      double m2Res   = mRes*mRes;
      double GamMRat = wRes/mRes;
      double sHat = phaseSpacePtr->sHat();
      double weightBW = pow2(sHat - m2Res) + pow2(sHat * GamMRat);      
      double m = std::sqrt(sHat)/((double)m_energyMode(settingsPtr)*1000.0);
 
      switch(m_energyMode(settingsPtr)){
      case 8:
    if(m < 0.0375){
      weightBW *= 121.88e-12*std::exp(13.0*m);
    }else{
      weightBW *= 1.0e-12*std::exp(18.5*m-1.4*std::log(m));
    }
 
    break;
 
      case 13:
    if(m < 0.023){
      weightBW *= 102.77e-12*std::exp(11.5*m);
    }else if(m < 0.231){
      weightBW *= 1.0e-12*std::exp(16.1*m-1.2*std::log(m));
    }else{
      weightBW *= 1.8675e-16*std::exp(31.7*m-4.6*std::log(m));
    }
    break;
 
      default:
    throw std::runtime_error("Unknown WprimeFlat:EnergyMode - should be either 8 or 13!");
 
      }
 
      return weightBW;
    }
 
  private:
    
    Pythia8_UserHooks::UserSetting<int> m_energyMode;
 
  };  
 
} // end namespace Pythia8