Commit a3d96688 authored by Frank Siegert's avatar Frank Siegert Committed by Ewelina Maria Lobodzinska
Browse files

Sherpa_i: Add base fragments for HF fusing.

parent a6691454
# taken from https://gitlab.com/sherpa-team/sherpa/merge_requests/228
# CAUTION: when including this, you have to set the following process-dependent setting in your JO:
# If two strong couplings are involved at Born level, `FUSING_DIRECT_FACTOR=1` (e.g. Zbb).
# If there are four such couplings, `FUSING_DIRECT_FACTOR=2` (e.g. ttbb).
genSeq.Sherpa_i.Parameters += [ "USERHOOK=Fusing_Direct",
"CSS_SCALE_SCHEME=2",
"CSS_EVOLUTION_SCHEME=3" ]
genSeq.Sherpa_i.PluginCode += r"""
#include "SHERPA/Tools/Userhook_Base.H"
#include "ATOOLS/Org/Message.H"
#include "SHERPA/Single_Events/Event_Handler.H"
#include "ATOOLS/Org/Data_Reader.H"
#include "ATOOLS/Phys/Weight_Info.H"
#include <algorithm>
#include "MODEL/Main/Running_AlphaS.H"
#include "SHERPA/Main/Sherpa.H"
#include "SHERPA/Initialization/Initialization_Handler.H"
#include "ATOOLS/Org/Data_Reader.H"
#include "CSSHOWER++/Main/CS_Shower.H"
#include "PDF/Main/Shower_Base.H"
#include "PHASIC++/Process/MCatNLO_Process.H"
/*
* counter-terms applied to the X+bb events in a fused sample.
*/
using namespace ATOOLS;
using namespace SHERPA;
class Fusing_Direct_Hook : public Userhook_Base {
private:
MODEL::Running_AlphaS *p_as;
Blob_List * p_bl;
Sherpa* p_sherpa;
double m_factor;
public:
Fusing_Direct_Hook(const Userhook_Arguments args) :
Userhook_Base("Example"), p_as(MODEL::as), p_sherpa(args.p_sherpa)
{
msg_Debugging()<<"Fusing_Direct Hook active."<<std::endl;
ATOOLS::Data_Reader *reader = p_sherpa->GetInitHandler()->DataReader();
m_factor = reader->GetValue<double>("FUSING_DIRECT_FACTOR", 2.);
}
~Fusing_Direct_Hook() {}
ATOOLS::Return_Value::code Run(ATOOLS::Blob_List* blobs, double &weight) {
p_bl = blobs;
String_BlobDataBase_Map bdmap = p_bl->FindFirst(btp::Signal_Process)->GetData();
auto search_mewinfo = bdmap.find("MEWeightInfo");
if (search_mewinfo==bdmap.end()) {
THROW(fatal_error,"No MEWeightinfo found in singnal blob!");
}
auto search_weight = bdmap.find("Weight");
if (search_weight==bdmap.end()) {
THROW(fatal_error,"No Weight found in singnal blob!");
}
ME_Weight_Info * me_w_info = search_mewinfo->second->Get<ME_Weight_Info *>();
double weight_bl = search_weight->second->Get<double>();
double mur2 = me_w_info->m_mur2;
ATOOLS::Flavour bquark = ATOOLS::Flavour(5);
double mb2 = bquark.Mass()*bquark.Mass();
double alphas((*p_as)(mur2));
double TR = 0.5;
double born_weight = me_w_info->m_B;
double new_weight(weight_bl);
double correction(0);
double sum_meweight = me_w_info->m_B + me_w_info->m_K + me_w_info->m_KP + me_w_info->m_VI;
//apply only to S-Events
if ( me_w_info->m_type ==(ATOOLS::mewgttype::code::METS |ATOOLS::mewgttype::code::H) ){
msg_Debugging() << "H-Event, skip alpha_s correction." << std::endl;
return Return_Value::Nothing;
}
PHASIC::Process_Vector procs = p_sherpa->GetInitHandler()->GetMatrixElementHandler()->AllProcesses();
if (procs.size()!=1)
THROW(fatal_error,"Too many procs!");
PHASIC::MCatNLO_Process * mcproc = dynamic_cast<PHASIC::MCatNLO_Process* >(procs.at(0));
if (mcproc==NULL){
THROW(fatal_error,"no MC@MLO process found! For use with separate LO-Process, use K-Factor!");
}
Cluster_Amplitude * ampl = mcproc->GetAmplitude();
double muf2 = ampl->KT2();
// ********** gg initial state ****************
if ((me_w_info->m_fl1 ==21) && (me_w_info->m_fl2 ==21)){
double l = log(mur2/muf2);
correction = alphas * 2.*TR/(3.*M_PI) * l ;
msg_Debugging() << "gg initial state. correction(" << correction <<")." << std::endl;
}
// ********** qq initial state ****************
if ((me_w_info->m_fl1 !=21) && (me_w_info->m_fl2 !=21)){
double l = log(mur2/mb2);
correction = alphas * 2.*TR/(3.*M_PI) * l ;
msg_Debugging() << "qq initial state. correction(" << correction <<")." << std::endl;
}
correction *= m_factor;
new_weight = weight_bl* (1. - correction* born_weight/sum_meweight);
weight = new_weight;
p_bl->FindFirst(btp::Signal_Process)->AddData("Weight",new Blob_Data<double>(new_weight));
// TODO: calculate counter-terms based on the muR variations. not done yet, since numerical impact is small.
auto search_varweights = bdmap.find("Variation_Weights");
if (search_varweights==bdmap.end()) {
THROW(fatal_error,"No VarWeight found in singnal blob!");
}
Variation_Weights var_weights= search_varweights->second->Get<Variation_Weights >();
var_weights*=(1. - correction* born_weight/sum_meweight);
p_bl->FindFirst(btp::Signal_Process)->AddData("Variation_Weights",new Blob_Data<Variation_Weights>(var_weights));
return Return_Value::Nothing;
}
void Finish() {
}
};
DECLARE_GETTER(Fusing_Direct_Hook,"Fusing_Direct",
Userhook_Base,Userhook_Arguments);
Userhook_Base *ATOOLS::Getter<Userhook_Base,Userhook_Arguments,Fusing_Direct_Hook>::
operator()(const Userhook_Arguments &args) const
{
return new Fusing_Direct_Hook(args);
}
void ATOOLS::Getter<Userhook_Base,Userhook_Arguments,Fusing_Direct_Hook>::
PrintInfo(std::ostream &str,const size_t width) const
{
str<<"Fusing_Direct userhook";
}
"""
# taken from https://gitlab.com/sherpa-team/sherpa/merge_requests/228
genSeq.Sherpa_i.Parameters += [ "USERHOOK=Fusing_Fragmentation",
"CSS_SCALE_SCHEME=2",
"CSS_EVOLUTION_SCHEME=3"
"FUSING_FRAGMENTATION_STORE_AS_WEIGHT=1" ]
genSeq.Sherpa_i.PluginCode += r"""
#include "SHERPA/Tools/Userhook_Base.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Math/Random.H"
#include "ATOOLS/Org/Run_Parameter.H"
#include "SHERPA/Single_Events/Event_Handler.H"
#include "ATOOLS/Org/Data_Reader.H"
#include "ATOOLS/Phys/Weight_Info.H"
#include <algorithm>
#include "MODEL/Main/Running_AlphaS.H"
#include "SHERPA/Main/Sherpa.H"
#include "SHERPA/Initialization/Initialization_Handler.H"
#include "ATOOLS/Org/Data_Reader.H"
/*
* implements the fragmentation component of the fusing algorithm based
* on a combined cluster amplitude.
*/
using namespace ATOOLS;
using namespace SHERPA;
class Fusing_Fragmentation_Hook : public Userhook_Base {
private:
MODEL::Running_AlphaS *p_as;
Blob_List * p_bl;
Sherpa* p_sherpa;
size_t m_tops_in_proc;
bool m_nlo, m_check_ho;
bool m_store;
public:
Fusing_Fragmentation_Hook(const Userhook_Arguments args) :
Userhook_Base("Fusing_Fragmentation"), p_as(MODEL::as), p_sherpa(args.p_sherpa)
{
msg_Debugging()<<"Fusing_Fragmentation Hook active."<<std::endl;
ATOOLS::Data_Reader *reader = p_sherpa->GetInitHandler()->DataReader();
m_nlo = reader->GetValue<int>("FUSING_FRAGMENTATION_NLO_Mode", 1);
m_store = reader->GetValue<int>("FUSING_FRAGMENTATION_STORE_AS_WEIGHT", 0);
m_check_ho = reader->GetValue<int>("FUSING_FRAGMENTATION_CHECK_HO", 1);
msg_Debugging()<<METHOD<< "NLO(" << m_nlo << ") \n }" << std::endl;
// count number of tops for ttbb in case of rare g->tt cluster configurations
m_tops_in_proc=0;
PHASIC::Process_Base * proc = p_sherpa->GetInitHandler()->GetMatrixElementHandler()->AllProcesses().at(0);
ATOOLS::Flavour_Vector flavs = proc->Flavours();
for (int i=2; i< flavs.size(); i++){
ATOOLS::Flavour flav = flavs.at(i);
if (flav.Kfcode()==6) m_tops_in_proc ++;
}
}
~Fusing_Fragmentation_Hook() {}
ATOOLS::Return_Value::code Run(ATOOLS::Blob_List* blobs, double &weight) {
p_bl = blobs;
ATOOLS::String_BlobDataBase_Map bdmap = p_bl->FindFirst(ATOOLS::btp::Shower)->GetData();
auto search = bdmap.find("AllAmplitudes");
if (search==bdmap.end()) {
THROW(fatal_error,"No matching amplitude found in blob. This algorithm works only with rel-2-2-7 or later!");
}
ATOOLS::Cluster_Amplitude * ampl = search->second->Get<ATOOLS::Cluster_Amplitude*>();
bool veto = true;
//find first amplitude with a bb-pair
size_t pos=0;
while (true){
if (ampl==NULL) break;
pos = CheckBPair(ampl);
if (pos==0){
ampl = ampl->Next();
} else break;
}
if (pos==0) veto=false; // no bb-pair found
else veto = !CheckFinalDir(ampl);
// redo veto for certain configurations (see code below)
if (veto && m_check_ho) {
if (CheckHigherOrder(ampl)){
msg_Debugging() << "skip Fusing Veto, configuration is of higher order.";
veto=false;
}
}
return ReturnFunc(!veto, weight);
}
Return_Value::code ReturnFunc(bool retval, double & event_weight){
//retval=false: veto Event
//retval=true: keep Event
String_BlobDataBase_Map bdmap = p_bl->FindFirst(btp::Signal_Process)->GetData();
auto search_weight = bdmap.find("Weight");
if (search_weight==bdmap.end()) {
THROW(fatal_error,"No Weight found in singnal blob!");
}
double weight_bl = search_weight->second->Get<double>();
if (retval==false){
msg_Debugging() << "false, veto configuration. \n";
}
// if not storing weights: do either Veto or not
if (!m_store) {
if (retval) return Return_Value::Nothing;
else return Return_Value::New_Event;
}
// if m_store: store additional weight(zero) in bloblist and slip veto
double new_weight = weight_bl;
if (retval==false) new_weight=0;
p_bl->FindFirst(ATOOLS::btp::Signal_Process)->AddData("UserHook",new ATOOLS::Blob_Data<double>(new_weight));
return Return_Value::Nothing;
}
/////////////////////// helper functions Veto
bool CheckFinalDir(ATOOLS::Cluster_Amplitude *ampl){
/*return only true, if a sufficient high number of light partons is present in this amplitude's final state
if the tops are clustered to a gluon before, there is one more light parton allowed in this step */
// true: keep amplitude
// false: do veto
size_t num_q(0);
size_t num_g(0);
size_t num_tops(0);
for (int i=2; i<ampl->Legs().size(); i++){
ATOOLS::Cluster_Leg *leg = ampl->Legs().at(i);
if (leg->Flav().IsGluon() && !leg->FromDec()) num_g++;
if (leg->Flav().Kfcode()<5 && !leg->FromDec()) num_q++;
if (leg->Flav().Kfcode()==6 && !leg->FromDec()) num_tops++;
}
if (m_nlo){
if (num_tops < m_tops_in_proc){
if ((num_q + num_g)>2) return true;
else return false;
}
else{
if ( (num_q + num_g)>1) return true;
else return false;
}
}
else{
if (num_tops < m_tops_in_proc){
if ((num_q >0) || (num_g>1)) return true;
else return false;
}
else{
if ((num_q >0) || (num_g>0)) return true;
else return false;
}
}
}
size_t CheckBPair(ATOOLS::Cluster_Amplitude *ampl){
if (ampl==NULL) return 0;
bool bfound(false), bbfound(false);
for (int i=2; i<ampl->Legs().size(); i++){
ATOOLS::Cluster_Leg *leg = ampl->Legs().at(i);
if ( !leg->FromDec() && leg->Flav().Kfcode()==5){
if (leg->Flav().IsAnti()) {
bbfound=true;
}else{
bfound=true;
}
}
}
if (bfound && bbfound) return 1;
return 0;
}
bool CheckHigherOrder (ATOOLS::Cluster_Amplitude *ampl){
/*check for configurations as bb->ttbb or bb->bb which should not be vetoed,
* since this is not part of 4F ttbb.
implemantation: redo veto, if bb-> x bb (j) and in x is no particle of the 93 container.
returns true, if such a config is found
*/
if (!CheckBPair(ampl)) {
THROW(fatal_error,"amplitude is missing b-quarks!");
}
msg_Debugging() << "Fusing, check higher order: " << ampl;
for (int i =0; i < 2; i++ ){
ATOOLS::Flavour flav = ampl->Leg(i)->Flav();
if (flav.Kfcode()!=5) return false;
}
return true;
}
void Finish() {
}
};
DECLARE_GETTER(Fusing_Fragmentation_Hook,"Fusing_Fragmentation",
Userhook_Base,Userhook_Arguments);
Userhook_Base *ATOOLS::Getter<Userhook_Base,Userhook_Arguments,Fusing_Fragmentation_Hook>::
operator()(const Userhook_Arguments &args) const
{
return new Fusing_Fragmentation_Hook(args);
}
void ATOOLS::Getter<Userhook_Base,Userhook_Arguments,Fusing_Fragmentation_Hook>::
PrintInfo(std::ostream &str,const size_t width) const
{
str<<"Fusing_Fragmentation_UserHook";
}
"""
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment