Ticket #317: OpenLoopsReplace_Interface.C

#include "PHASIC++/Process/Process_Base.H"
#include "PHASIC++/Process/ME_Generator_Base.H"
#include <dirent.h>


#include "PHASIC++/Process/Tree_ME2_Base.H"
#include "PHASIC++/Process/Virtual_ME2_Base.H"

namespace MODEL {
  class Model_Base;
}

namespace OpenLoopsReplace {

typedef void (*Amp2Func)(double* moms, double* M2L0, double* M2L1, double* IRL1,
                         double* M2L2, double* IRL2);

typedef void (*PermutationFunc)(int* permutation);




  class OpenLoopsReplace_Particle {
  public:
    std::string m_name, m_faname;
    size_t m_order;
    OpenLoopsReplace_Particle() : m_name(""), m_faname(""), m_order(0) {}
    OpenLoopsReplace_Particle(std::string name, std::string faname, size_t order) :
      m_name(name), m_faname(faname), m_order(order) {}
  };


  class OpenLoopsReplace_Interface: public PHASIC::ME_Generator_Base {

    void InitializeParticles();

    static std::string MatchOptions(std::vector<std::string> options,
                                    int oew, int oqcd, int nloop);
    static int SelectInfo(const dirent *entry);
    static bool Order(const std::pair<size_t, ATOOLS::Flavour> &a,
                      const std::pair<size_t, ATOOLS::Flavour> &b);

  public :
    OpenLoopsReplace_Interface() : ME_Generator_Base("OpenLoopsReplace") {}
    ~OpenLoopsReplace_Interface() {}

    bool Initialize(const std::string &path,const std::string &file,
                    MODEL::Model_Base *const model,
                    BEAM::Beam_Spectra_Handler *const beam,
                    PDF::ISR_Handler *const isr);

    static void SetParameter(const std::string & key, double value);
    static void SetParameter(const std::string & key, int value);
    static void SetParameter(const std::string & key, std::string value);
    static double GetDoubleParameter(const std::string & key);
    static int GetIntParameter(const std::string & key);
    static void FlushParameters();

    static std::string GetProcessPermutation(const ATOOLS::Flavour_Vector& flavs_orig,
                                             std::vector<int>& permutation);
    static std::pair<std::string, std::string> ScanFiles(std::string& process,
                                                         int oew, int oqcd, int nloop);
    static ATOOLS::Flavour_Vector MapFlavours(const ATOOLS::Flavour_Vector& orig);

    PHASIC::Process_Base *InitializeProcess(const PHASIC::Process_Info &pi, bool add)
    { return NULL; }
    int  PerformTests() { return 1; }
    bool NewLibraries() { return false; }
    void SetClusterDefinitions(PDF::Cluster_Definitions_Base *const defs) {}
    ATOOLS::Cluster_Amplitude *ClusterConfiguration
    (PHASIC::Process_Base *const proc,const size_t &mode) { return NULL; }

    static std::string s_olprefix;
    static std::vector<std::string> s_allowed_libs;
    static bool s_generate_list;
    static std::set<std::string> s_shoppinglist;
    static std::map<int, OpenLoopsReplace_Particle> s_particles;
  };


  class OpenLoopsReplace_Born: public PHASIC::Tree_ME2_Base {

    Amp2Func m_amp2;
    PermutationFunc m_permutationfunc;
    std::vector<int> m_permutation;
    static OpenLoopsReplace_Interface* s_interface;

    size_t m_oqcd, m_oew;
    double m_symfac;

  public:

    OpenLoopsReplace_Born(const PHASIC::Process_Info& pi,
                   const ATOOLS::Flavour_Vector& flavs,
                   Amp2Func amp2,
                   PermutationFunc permutationfunc,
                   std::vector<int> permutation,
                   std::string functag);

    ~OpenLoopsReplace_Born();

    double Calc(const ATOOLS::Vec4D_Vector& momenta);

    inline static void SetInterface(OpenLoopsReplace_Interface* interface)
    {
      s_interface=interface;
    }
  };


  class OpenLoopsReplace_Virtual : public PHASIC::Virtual_ME2_Base {

    Amp2Func m_amp2;
    PermutationFunc m_permutationfunc;
    std::vector<int> m_permutation;
    static OpenLoopsReplace_Interface* s_interface;
    int m_oew, m_oqcd;
    int m_write_points;
    std::ofstream m_points_file;

    bool     m_analyse;
    std::map<std::string,ATOOLS::Histogram *> m_histograms;

  public:

    OpenLoopsReplace_Virtual(const PHASIC::Process_Info& pi,
                      const ATOOLS::Flavour_Vector& flavs,
                      Amp2Func amp2,
                      PermutationFunc permutationfunc,
                      std::vector<int> permutation,
                      std::string functag);

    ~OpenLoopsReplace_Virtual();
  
    inline bool SetColours(const ATOOLS::Vec4D_Vector& momenta) {
      return true;
    }

    inline double Eps_Scheme_Factor(const ATOOLS::Vec4D_Vector& mom) {
      return 4.*M_PI;
    }
  
    void Calc(const ATOOLS::Vec4D_Vector& momenta);

    inline static void SetInterface(OpenLoopsReplace_Interface* interface)
    {
      s_interface=interface;
    }

  };

}



#include "ATOOLS/Org/CXXFLAGS_PACKAGES.H"
#include "MODEL/Main/Model_Base.H"
#include "PHASIC++/Main/Phase_Space_Handler.H"
#include "ATOOLS/Org/Data_Reader.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Org/MyStrStream.H"
#include "ATOOLS/Org/Library_Loader.H"
#include "ATOOLS/Org/Run_Parameter.H"
#include <algorithm>
#include <sys/stat.h>

using namespace PHASIC;
using namespace MODEL;
using namespace ATOOLS;
using namespace std;


extern "C" {

  void openloops_welcome_(char* str, int* len_str);

  void loop_parameters_write_();

  void ol_getparameter_double_c_(const char* key, double* val, int* err);
  void ol_getparameter_int_c_(const char* key, int* val, int* err);
  void ol_setparameter_double_c_(const char* key, double* val, int* err);
  void ol_setparameter_int_c_(const char* key, int* val, int* err);
  void ol_setparameter_string_c_(const char* key, const char* val, int* err);
  void ol_parameters_flush_();

}


namespace OpenLoopsReplace {

  std::map<int, OpenLoopsReplace_Particle> OpenLoopsReplace_Interface::s_particles;
  std::string OpenLoopsReplace_Interface::s_olprefix;
  std::vector<std::string> OpenLoopsReplace_Interface::s_allowed_libs;
  bool OpenLoopsReplace_Interface::s_generate_list;
  set<string> OpenLoopsReplace_Interface::s_shoppinglist;

  bool OpenLoopsReplace_Interface::Initialize(const string &path,const string &file,
                                       MODEL::Model_Base *const model,
                                       BEAM::Beam_Spectra_Handler *const beam,
                                       PDF::ISR_Handler *const isr)
  {
    InitializeParticles();

    // find OL installation prefix with several overwrite options
    struct stat st;
    Data_Reader reader(" ",";","#","=");
    s_olprefix = rpa->gen.Variable("SHERPA_CPP_PATH")+"/Process/OpenLoopsReplace";
    if(stat(s_olprefix.c_str(),&st) != 0) s_olprefix = OPENLOOPS_PREFIX;
    s_olprefix = reader.GetValue<string>("OL_PREFIX", s_olprefix);
    msg_Info()<<"Initialising OpenLoopsReplace generator from "<<s_olprefix<<endl;

    // load library dynamically
    s_loader->AddPath(s_olprefix+"/lib");
    s_loader->AddPath(s_olprefix+"/proclib");
    if (!s_loader->LoadLibrary("openloops")) THROW(fatal_error, "Failed to load libopenloops.");

    // set particle masses/widths
    int tmparr[] = {kf_e, kf_mu, kf_tau, kf_u, kf_d, kf_s, kf_c, kf_b, kf_t, kf_Wplus, kf_Z, kf_h0};
    vector<int> pdgids (tmparr, tmparr + sizeof(tmparr) / sizeof(tmparr[0]) );
    int dummy;
    for (size_t i=0; i<pdgids.size(); ++i) {
      if (Flavour(pdgids[i]).Mass()>0.0) SetParameter("mass("+ToString(pdgids[i])+")", Flavour(pdgids[i]).Mass());
      if (Flavour(pdgids[i]).Width()>0.0) SetParameter("width("+ToString(pdgids[i])+")", Flavour(pdgids[i]).Width());
    }

    // set remaining OL parameters specified by user
    vector<string> parameters;
    reader.VectorFromFile(parameters,"OL_PARAMETERS");
    for (size_t i=1; i<parameters.size(); i=i+2) SetParameter(parameters[i-1], parameters[i]);
    ol_parameters_flush_();

    reader.VectorFromFile(s_allowed_libs, "OL_ALLOWED_LIBS");
    s_generate_list=reader.GetValue<size_t>("OL_GENERATE_LIST", false);
    if (s_allowed_libs.size()) PRINT_VAR(s_allowed_libs);

    int lenws=1200;
    char welcomestr[lenws];
    openloops_welcome_(welcomestr, &lenws);
    msg_Info()<<std::string(welcomestr,lenws)<<std::endl;

    MyStrStream cite;
    cite<<"The OpenLoopsReplace library~\\cite{Cascioli:2011va} of virtual"<<endl
        <<"matrix elements has been used. "<<endl;
    if (GetIntParameter("redlib1")==1 || GetIntParameter("redlib1")==7 ||
        GetIntParameter("redlib2")==1 || GetIntParameter("redlib2")==7) {
      cite<<"It is partly based on the tensor integral reduction described "<<endl
          <<"in~\\cite{Denner:2002ii,Denner:2005nn,Denner:2010tr}."<<endl;
    }
    if (GetIntParameter("redlib1")==5 || GetIntParameter("redlib2")==5) {
      cite<<"It is partly based on the integrand reduction described "<<endl
          <<"in~\\cite{Ossola:2007ax,vanHameren:2010cp}."<<endl;
    }
    if (GetIntParameter("redlib1")==6 || GetIntParameter("redlib2")==6) {
      cite<<"It is partly based on the integrand reduction described "<<endl
          <<"in~\\cite{Mastrolia:2010nb,vanHameren:2010cp}."<<endl;
    }
    if (GetIntParameter("redlib1")==8 || GetIntParameter("redlib2")==8) {
      cite<<"It is partly based on the integrand reduction Ninja."<<endl;
    }
    rpa->gen.AddCitation(1,cite.str());

    OpenLoopsReplace_Virtual::SetInterface(this);
    OpenLoopsReplace_Born::SetInterface(this);
    return true;
  }

  void OpenLoopsReplace_Interface::InitializeParticles()
  {
    s_particles.clear();
    s_particles[kf_nue] = OpenLoopsReplace_Particle("ne", "F[1,{1}]", 0);
    s_particles[-kf_nue] = OpenLoopsReplace_Particle("nex", "-F[1,{1}]", 1);
    s_particles[kf_numu] = OpenLoopsReplace_Particle("nm", "F[1,{2}]", 2);
    s_particles[-kf_numu] = OpenLoopsReplace_Particle("nmx", "-F[1,{2}]", 3);
    s_particles[kf_nutau] = OpenLoopsReplace_Particle("nl", "F[1,{3}]", 4);
    s_particles[-kf_nutau] = OpenLoopsReplace_Particle("nlx", "-F[1,{3}]", 5);
    s_particles[kf_e] = OpenLoopsReplace_Particle("e", "F[2,{1}]", 6);
    s_particles[-kf_e] = OpenLoopsReplace_Particle("ex", "-F[2,{1}]", 7);
    s_particles[kf_mu] = OpenLoopsReplace_Particle("m", "F[2,{2}]", 8);
    s_particles[-kf_mu] = OpenLoopsReplace_Particle("mx", "-F[2,{2}]", 9);
    s_particles[kf_tau] = OpenLoopsReplace_Particle("l", "F[2,{3}]", 10);
    s_particles[-kf_tau] = OpenLoopsReplace_Particle("lx", "-F[2,{3}]", 11);
    s_particles[kf_u] = OpenLoopsReplace_Particle("u", "F[3,{1}]", 12);
    s_particles[-kf_u] = OpenLoopsReplace_Particle("ux", "-F[3,{1}]", 13);
    s_particles[kf_c] = OpenLoopsReplace_Particle("c", "F[3,{2}]", 14);
    s_particles[-kf_c] = OpenLoopsReplace_Particle("cx", "-F[3,{2}]", 15);
    s_particles[kf_t] = OpenLoopsReplace_Particle("t", "F[3,{3}]", 16);
    s_particles[-kf_t] = OpenLoopsReplace_Particle("tx", "-F[3,{3}]", 17);
    s_particles[kf_d] = OpenLoopsReplace_Particle("d", "F[4,{1}]", 18);
    s_particles[-kf_d] = OpenLoopsReplace_Particle("dx", "-F[4,{1}]", 19);
    s_particles[kf_s] = OpenLoopsReplace_Particle("s", "F[4,{2}]", 20);
    s_particles[-kf_s] = OpenLoopsReplace_Particle("sx", "-F[4,{2}]", 21);
    s_particles[kf_b] = OpenLoopsReplace_Particle("b", "F[4,{3}]", 22);
    s_particles[-kf_b] = OpenLoopsReplace_Particle("bx", "-F[4,{3}]", 23);
    s_particles[kf_h0] = OpenLoopsReplace_Particle("h", "S[1]", 24);
    s_particles[kf_photon] = OpenLoopsReplace_Particle("a", "V[1]", 25);
    s_particles[kf_Z] = OpenLoopsReplace_Particle("z", "V[2]", 26);
    s_particles[-kf_Wplus] = OpenLoopsReplace_Particle("w", "V[3]", 27);
    s_particles[kf_Wplus] = OpenLoopsReplace_Particle("wx", "-V[3]", 28);
    s_particles[kf_gluon] = OpenLoopsReplace_Particle("g", "V[5]", 29);
  }

  std::string OpenLoopsReplace_Interface::MatchOptions(vector<string> options, int oew, int oqcd, int nloop) {
    for (size_t i=2; i<options.size(); ++i) {
      string option=options[i].substr(0, options[i].find("="));
      string value=options[i].substr(options[i].find("=")+1);

      if (option=="EW" && value!=ToString(oew)+",0") return "0";
      if (option=="QCD" && value!=ToString(oqcd)+","+ToString(nloop)) return "0";
      if (option=="CKMORDER") {
        int ckmorder=ToType<int>(value);
        if (ckmorder<3) {
          if (s_model->ComplexMatrixElement("CKM", 0,2)!=Complex(0.0,0.0) ||
              s_model->ComplexMatrixElement("CKM", 2,0)!=Complex(0.0,0.0)) {
            return "0";
          }
        }
        if (ckmorder<2) {
          if (s_model->ComplexMatrixElement("CKM", 1,2)!=Complex(0.0,0.0) ||
              s_model->ComplexMatrixElement("CKM", 2,1)!=Complex(0.0,0.0)) {
            return "0";
          }
        }
        if (ckmorder<1) {
          if (s_model->ComplexMatrixElement("CKM", 0,1)!=Complex(0.0,0.0) ||
              s_model->ComplexMatrixElement("CKM", 1,0)!=Complex(0.0,0.0)) {
            return "0";
          }
        }
      }
      const string nfs("nf");
      double nfsv=GetIntParameter(nfs);
      if (option=="nf" && ToType<int>(value)!=nfsv) return "0";
      if (option=="MD" && (double)GetDoubleParameter("mass(1)")>0.0) return "0";
      if (option=="MU" && (double)GetDoubleParameter("mass(2)")>0.0) return "0";
      if (option=="MS" && (double)GetDoubleParameter("mass(3)")>0.0) return "0";
      if (option=="MC" && (double)GetDoubleParameter("mass(4)")>0.0) return "0";
      if (option=="MB" && (double)GetDoubleParameter("mass(5)")>0.0) return "0";
      if (option=="MT" && (double)GetDoubleParameter("mass(6)")>0.0) return "0";
      if (option=="ME" && (double)GetDoubleParameter("mass(11)")>0.0) return "0";
      if (option=="MM" && (double)GetDoubleParameter("mass(13)")>0.0) return "0";
      if (option=="MT" && (double)GetDoubleParameter("mass(15)")>0.0) return "0";

      if (option=="map") return value;
    }
    return "1";
  }


  int OpenLoopsReplace_Interface::SelectInfo(const dirent *entry)
  {
    string fname(entry->d_name);
    if (fname.find(".info")!=string::npos) {
      return true;
    }
    else return false;
  }

  pair<string, string> OpenLoopsReplace_Interface::ScanFiles(string& process, int oew, int oqcd, int nloop)
  {
    struct dirent **entries;
    string procdatapath=s_olprefix+"/proclib";
    int n(scandir(procdatapath.c_str(),&entries,&SelectInfo,alphasort));
    if (n<0) THROW(fatal_error, "OpenLoopsReplace process dir "+procdatapath+" not found.");
    vector<string> files;
    for (int ifile=0; ifile<n; ++ifile) {
      files.push_back(string(entries[ifile]->d_name));
      free(entries[ifile]);
    }
    free(entries);

    for (int ifile=0; ifile<files.size(); ++ifile) {
      Data_Reader reader(" ",";","#","");
      reader.SetAddCommandLine(false);
      reader.SetInputFile(procdatapath+"/"+files[ifile]);
      reader.SetMatrixType(mtc::transposed);
      vector<vector<string> > content;
      reader.MatrixFromFile(content);
      for (size_t i=0; i<content.size(); ++i) {

        bool autoload=true;
        if (content[i].size()>0 && content[i][0]=="options") {
          for (size_t j=0; j<content[i].size(); ++j) {
            if (content[i][j]=="autoload=0") {
              autoload=false;
            }
          }
        }

        if (content[i].size()>2 && content[i][1]==process) {
          string match=MatchOptions(content[i], oew, oqcd, nloop);
          if (match=="0") {
            PRINT_INFO("Ignoring process with incompatible options.");
            continue;
          }
          else if (match!="1") {
            PRINT_INFO("Mapping "<<process<<" to "<<match<<".");
            process=match;
            return ScanFiles(process, oew, oqcd, nloop);
          }
          string grouptag=content[i][0];
          string process_subid=content[i][2];
          if (s_allowed_libs.size()>0 || autoload==false) {
            bool allowed=false;
            for (size_t i=0; i<s_allowed_libs.size(); ++i) {
              if (grouptag==s_allowed_libs[i]) allowed=true;
            }
            if (!allowed) continue;
          }
          return make_pair(grouptag, process_subid);
        }
      }
    }
    PRINT_INFO("Didn't find info file matching process "<<process);
    return make_pair("", "0");
  }


  bool OpenLoopsReplace_Interface::Order(const pair<size_t, Flavour> &a,const pair<size_t, Flavour> &b)
  {
    if (s_particles.find(a.second.HepEvt())==s_particles.end() ||
        s_particles.find(b.second.HepEvt())==s_particles.end()) {
      THROW(fatal_error, "Unknown particle.");
    }
    return s_particles[a.second.HepEvt()].m_order<s_particles[b.second.HepEvt()].m_order;
  }


  string OpenLoopsReplace_Interface::GetProcessPermutation(const Flavour_Vector& flavs_orig,
                                                    vector<int>& permutation) {
    vector<pair<size_t, Flavour> > flavs_ol(flavs_orig.size());
    for (size_t i=0; i<2; ++i) {
      flavs_ol[i]=make_pair(i,flavs_orig[i]);
    }
    for (size_t i=2; i<flavs_orig.size(); ++i) {
      flavs_ol[i]=make_pair(i,flavs_orig[i].Bar());
    }

    sort(flavs_ol.begin(), flavs_ol.end(), Order);

    permutation.clear();
    permutation.resize(flavs_ol.size());
    for (size_t i=0; i<flavs_ol.size(); ++i) {
      for (size_t j=0; j<flavs_ol.size(); ++j) {
        if (i==flavs_ol[j].first) {
          permutation[i]=int(j+1); // +1 for fortran openloops
          break;
        }
      }
    }

    string process;
    for (size_t i=0; i<flavs_ol.size(); ++i) {
      process += s_particles[flavs_ol[i].second.HepEvt()].m_name;
    }
    return process;
  }

  bool SortByFirstDec(pair<size_t, size_t> p1, pair<size_t, size_t> p2) {
    return p1.first>p2.first;
  }

  Flavour_Vector OpenLoopsReplace_Interface::MapFlavours(const Flavour_Vector& orig)
  {
    // for (size_t i=2; i<orig.size(); ++i) {
    //   if (false && orig[i].Width()!=0.0) {
    //     THROW(fatal_error, "Non-zero width of final state particle.");
    //   }
    // }

    /* Concept:
      For each family i=0,...,2:
      (1) Given a final state, determine the four (anti)lepton/neutrino
          multiplicities in the given process:
          a_nubar, a_nu, a_lbar, a_l

      (2) Compute the discriminant N[i] as
          N[i] = Ngen - (i+1) + Ngen*(a_nubar + a_nu*Nmax + a_lbar*Nmax^2 + a_l*Nmax^3)
          where Nmax should be chosen such that a_...<=Nmax.
          In practice one can safely set Nmax=10 and it will work for any
          process with <= 20 final-state leptons.
          It is also convenient to set Ngen=10, although i runs only from 1 to 3.

      (3) Reassign the lepton generations with a permutation
          p1 -> 1, p2 -> 2, p3 -> 3  such that  N[p1] > N[p2] > N[p3] */

    multiset<int> hepevt;
    for (size_t i=0; i<orig.size(); ++i) { hepevt.insert(orig[i].HepEvt()); }

    size_t Ngen(10), Nmax(10);
    vector<pair<size_t, size_t> > N(3);
    for (size_t i=0; i<3; ++i) {
      int nu_gen=10+2*(i+1);
      int l_gen=9+2*(i+1);

      size_t a_nu=hepevt.count(nu_gen);
      size_t a_nubar=hepevt.count(-nu_gen);
      size_t a_l=hepevt.count(l_gen);
      size_t a_lbar=hepevt.count(-l_gen);

      N[i]=make_pair(Ngen-(i+1) + Ngen*
                     (a_nubar+a_nu*Nmax+a_lbar*Nmax*Nmax+a_l*Nmax*Nmax*Nmax),
                     i);
    }

    sort(N.begin(), N.end(), SortByFirstDec);

    Flavour_Vector ret(orig);
    for (size_t i=0; i<3; ++i) {
      int nu_gen=10+2*(N[i].second+1);
      int nu_gen_new=10+2*(i+1);
      int l_gen=9+2*(N[i].second+1);
      int l_gen_new=9+2*(i+1);

      for (size_t j=0; j<orig.size(); ++j) {
        if (orig[j].Kfcode()==nu_gen) {
          ret[j]=Flavour(nu_gen_new, orig[j].IsAnti());
        }
        if (orig[j].Kfcode()==l_gen) {
          ret[j]=Flavour(l_gen_new, orig[j].IsAnti());
        }
      }
    }
    return ret;
  }

  double OpenLoopsReplace_Interface::GetDoubleParameter(const std::string & key) {
    int err(0);
    double value;
    ol_getparameter_double_c_(key.c_str(), &value, &err);
    return value;
  }
  int OpenLoopsReplace_Interface::GetIntParameter(const std::string & key) {
    int err(0);
    int value;
    ol_getparameter_int_c_(key.c_str(), &value, &err);
    return value;
  }
  template <class ValueType>
  void HandleParameterStatus(int err, const std::string & key, ValueType value) {
    if (err==0) {
      msg_Debugging()<<"Setting OpenLoopsReplace parameter: "<<key<<" = "<<value<<endl;
    }
    else if (err==1) {
      THROW(fatal_error, "Unknown OpenLoopsReplace parameter: "+key+" = "+ToString(value));
    }
    else if (err==2) {
      THROW(fatal_error, "Error setting OpenLoopsReplace parameter: "+key+" = "+ToString(value));
    }
  }
  void OpenLoopsReplace_Interface::SetParameter(const std::string & key, double value) {
    int err(0);
    ol_setparameter_double_c_(key.c_str(), &value, &err);
    HandleParameterStatus(err, key, value);
  }
  void OpenLoopsReplace_Interface::SetParameter(const std::string & key, int value) {
    int err(0);
    ol_setparameter_int_c_(key.c_str(), &value, &err);
    HandleParameterStatus(err, key, value);
  }
  void OpenLoopsReplace_Interface::SetParameter(const std::string & key, std::string value) {
    int err(0);
    ol_setparameter_string_c_(key.c_str(), value.c_str(), &err);
    HandleParameterStatus(err, key, value);
  }
  void OpenLoopsReplace_Interface::FlushParameters() {
    ol_parameters_flush_();
  }


}

using namespace OpenLoopsReplace;

DECLARE_GETTER(OpenLoopsReplace_Interface,"OpenLoopsReplace",ME_Generator_Base,ME_Generator_Key);

ME_Generator_Base *ATOOLS::Getter<ME_Generator_Base,ME_Generator_Key,
                                  OpenLoopsReplace_Interface>::
operator()(const ME_Generator_Key &key) const
{
  return new OpenLoopsReplace::OpenLoopsReplace_Interface();
}

void ATOOLS::Getter<ME_Generator_Base,ME_Generator_Key,OpenLoopsReplace_Interface>::
PrintInfo(ostream &str,const size_t width) const
{ 
  str<<"Interface to the OpenLoopsReplace loop ME generator"; 
}




#include "ATOOLS/Org/CXXFLAGS.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Org/Library_Loader.H"

using namespace PHASIC;
using namespace MODEL;
using namespace ATOOLS;
using namespace std;

namespace OpenLoopsReplace {

OpenLoopsReplace_Interface* OpenLoopsReplace_Born::s_interface=NULL;

OpenLoopsReplace_Born::OpenLoopsReplace_Born(const Process_Info& pi,
                               const Flavour_Vector& flavs,
                               Amp2Func amp2,
                               PermutationFunc permutationfunc,
                               std::vector<int> permutation,
                               std::string functag) :
  Tree_ME2_Base(pi, flavs),
  m_amp2(amp2), m_permutationfunc(permutationfunc),
  m_permutation(permutation)
{
  m_symfac=pi.m_fi.FSSymmetryFactor();
  m_symfac*=pi.m_ii.ISSymmetryFactor();
  m_oew=pi.m_oew;
  m_oqcd=pi.m_oqcd;
}

OpenLoopsReplace_Born::~OpenLoopsReplace_Born()
{
}

double OpenLoopsReplace_Born::Calc(const Vec4D_Vector& momenta)
{
  Vec4D_Vector m_moms(momenta);

  s_interface->SetParameter("alpha", AlphaQED());
  s_interface->SetParameter("alphas", AlphaQCD());
  s_interface->FlushParameters();

  double M2L0;
  std::vector<double> M2L1(3), M2L2(5), IRL1(3), IRL2(5);

  m_permutationfunc(&m_permutation[0]);
  m_amp2(&m_moms[0][0], &M2L0, &M2L1[0], &IRL1[0], &M2L2[0], &IRL2[0]);

  if (IsZero(M2L1[1]) && IsZero(M2L1[2]) &&
      IsZero(IRL1[1]) && IsZero(IRL1[2]) &&
      IsZero(M2L2[1]) && IsZero(M2L2[2]) && IsZero(M2L2[3]) && IsZero(M2L2[4]) &&
      IsZero(IRL2[1]) && IsZero(IRL2[2]) && IsZero(IRL2[3]) && IsZero(IRL2[4])) {
    // OL returns ME2 including 1/symfac, but Calc is supposed to return it
    // without 1/symfac, thus multiplying with symfac here
    if (IsZero(M2L0)) return m_symfac*M2L2[0];
    else return m_symfac*M2L0;
  }
  else {
    PRINT_INFO("Poles non-zero. Returning 0.");
    return 0.0;
  }
}

}

using namespace OpenLoopsReplace;

DECLARE_TREEME2_GETTER(OpenLoopsReplace_Born,"OpenLoopsReplace_Born")
Tree_ME2_Base *ATOOLS::Getter<Tree_ME2_Base,Process_Info,OpenLoopsReplace_Born>::
operator()(const Process_Info &pi) const
{
  DEBUG_FUNC(pi);
  if (pi.m_loopgenerator!="OpenLoopsReplace") return NULL;
  if (pi.m_fi.m_nloewtype!=nlo_type::lo) return NULL;
  if (pi.m_fi.m_nloqcdtype!=nlo_type::lo &&
      pi.m_fi.m_nloqcdtype!=nlo_type::born &&
      pi.m_fi.m_nloqcdtype!=nlo_type::real) return NULL;
  if (MODEL::s_model->Name()!="SM") return NULL;

  Flavour_Vector flavs=pi.ExtractFlavours();
  Flavour_Vector map_flavs=OpenLoopsReplace_Interface::MapFlavours(flavs);
  msg_Tracking()<<endl<<flavs<<" --> "<<map_flavs<<endl;

  vector<int> permutation;
  string process=OpenLoopsReplace_Interface::GetProcessPermutation(map_flavs, permutation);
  pair<string, string> groupsub=OpenLoopsReplace_Interface::ScanFiles(process, pi.m_oew, pi.m_oqcd, 0);
  string grouptag=groupsub.first;
  string subid=groupsub.second;
  if (grouptag!="") {
    // symbols in fortran are always defined as lower case
    string lc_functag(grouptag+"_"+process+"_"+subid+"_");
    for (size_t i(0);i<lc_functag.length();++i)
      lc_functag[i]=tolower(lc_functag[i]);
    vector<string> suffixes;
    suffixes.push_back("1sL"); // deprecated
    suffixes.push_back("ls");
    suffixes.push_back("lst");
    suffixes.push_back("lps");
    suffixes.push_back("lpst");
    void *ampfunc, *permfunc;
    for (size_t i=0; i<suffixes.size(); ++i) {
      string libraryfile="openloops_"+grouptag+"_"+suffixes[i];
      ampfunc=s_loader->GetLibraryFunction(libraryfile,"vamp2_"+lc_functag);
      permfunc=s_loader->GetLibraryFunction(libraryfile,"set_permutation_"+lc_functag);
      if (ampfunc!=NULL && permfunc!=NULL) break;
    }
    if (ampfunc==NULL || permfunc==NULL) {
      PRINT_INFO("Didn't find functions");
      return NULL;
    }

    msg_Info()<<endl;
    PRINT_INFO("Initialising OpenLoopsReplace Born for "<<flavs<<": "<<lc_functag);
    return new OpenLoopsReplace_Born(pi, flavs, (Amp2Func) ampfunc,
                              (PermutationFunc) permfunc, permutation, lc_functag);
  }
  else {
    return NULL;
  }
}


#include "MODEL/Main/Model_Base.H"
#include "ATOOLS/Org/Run_Parameter.H"
#include "ATOOLS/Org/Exception.H"
#include "ATOOLS/Org/MyStrStream.H"
#include "ATOOLS/Org/Data_Reader.H"
#include "ATOOLS/Org/Library_Loader.H"
#include "ATOOLS/Math/Poincare.H"


using namespace PHASIC;
using namespace ATOOLS;
using namespace std;

namespace OpenLoopsReplace {

OpenLoopsReplace_Interface* OpenLoopsReplace_Virtual::s_interface=NULL;

OpenLoopsReplace_Virtual::OpenLoopsReplace_Virtual(const Process_Info& pi,
                                     const Flavour_Vector& flavs,
                                     Amp2Func amp2,
                                     PermutationFunc permutationfunc,
                                     std::vector<int> permutation,
                                     std::string functag) :
  Virtual_ME2_Base(pi, flavs),
  m_amp2(amp2), m_permutationfunc(permutationfunc),
  m_permutation(permutation), m_analyse(true)
{
  m_oew=pi.m_oew;
  m_oqcd=pi.m_oqcd;
}

OpenLoopsReplace_Virtual::~OpenLoopsReplace_Virtual()
{
}


void OpenLoopsReplace_Virtual::Calc(const Vec4D_Vector& momenta) {
  Vec4D_Vector m_moms(momenta);

  s_interface->SetParameter("alpha", AlphaQED());
  s_interface->SetParameter("alphas", AlphaQCD());
  s_interface->SetParameter("mu", sqrt(m_mur2));
  s_interface->FlushParameters();
  
  double M2L0;
  vector<double> M2L1(3), M2L2(5), IRL1(3), IRL2(5);

  MyTiming* timing;
  if (msg_LevelIsDebugging()) {
    timing = new MyTiming();
    timing->Start();
  }
  m_permutationfunc(&m_permutation[0]);
  m_amp2(&m_moms[0][0], &M2L0, &M2L1[0], &IRL1[0], &M2L2[0], &IRL2[0]);
  if (msg_LevelIsDebugging()) {
    timing->Stop();
    PRINT_INFO(momenta[2][0]<<" "<<m_flavs<<" user="<<timing->UserTime()
               <<" real="<<timing->RealTime()<<" sys="<<timing->SystemTime());
  }

  double B(M2L0), V_finite(M2L1[0]), V_eps(M2L1[1]), V_eps2(M2L1[2]);

  // factor which by Sherpa convention has to be divided out at this stage
  double factor=B*AlphaQCD()/2.0/M_PI;

  m_born=B;
  m_res.Finite()=(V_finite/factor);
  m_res.IR()=(V_eps/factor);
  m_res.IR2()=(V_eps2/factor);
}

}


using namespace OpenLoopsReplace;

void dummyamp2func(double* moms, double* M2L0, double* M2L1, double* IRL1, double* M2L2, double* IRL2)
{ THROW(normal_exit, "Shopping list generated."); }
void dummypermfunc(int* permutation)
{ THROW(normal_exit, "Shopping list generated."); }

DECLARE_VIRTUALME2_GETTER(OpenLoopsReplace_Virtual,"OpenLoopsReplace_Virtual")
Virtual_ME2_Base *ATOOLS::Getter<Virtual_ME2_Base,Process_Info,OpenLoopsReplace_Virtual>::
operator()(const Process_Info &pi) const
{
  DEBUG_FUNC(pi);
  if (pi.m_loopgenerator!="OpenLoopsReplace") return NULL;
  if (pi.m_fi.m_nloewtype!=nlo_type::lo) return NULL;
  if (pi.m_fi.m_nloqcdtype!=nlo_type::loop) return NULL;
  if (MODEL::s_model->Name()!="SM") return NULL;

  Flavour_Vector flavs=pi.ExtractFlavours();
  Flavour_Vector map_flavs=OpenLoopsReplace_Interface::MapFlavours(flavs);
  msg_Tracking()<<endl<<flavs<<" --> "<<map_flavs<<endl;

  vector<int> permutation;
  string process=OpenLoopsReplace_Interface::GetProcessPermutation(map_flavs, permutation);
  pair<string, string> groupsub=OpenLoopsReplace_Interface::ScanFiles(process, pi.m_oew, pi.m_oqcd-1, 1);
  string grouptag=groupsub.first;
  string subid=groupsub.second;
  if (grouptag!="") {
    // symbols in fortran are always defined as lower case
    string lc_functag(grouptag+"_"+process+"_"+subid+"_");
    for (size_t i(0);i<lc_functag.length();++i)
      lc_functag[i]=tolower(lc_functag[i]);
    vector<string> suffixes;
    suffixes.push_back("1tL"); // deprecated
    suffixes.push_back("1L"); // deprecated
    suffixes.push_back("1ptL"); // deprecated
    suffixes.push_back("0L"); // deprecated
    suffixes.push_back("l");
    suffixes.push_back("lt");
    suffixes.push_back("lpt");
    suffixes.push_back("ls");
    suffixes.push_back("lp");
    suffixes.push_back("lst");
    suffixes.push_back("lps");
    suffixes.push_back("lpst");
    void *ampfunc, *permfunc;
    for (size_t i=0; i<suffixes.size(); ++i) {
      string libraryfile="openloops_"+grouptag+"_"+suffixes[i];
      ampfunc=s_loader->GetLibraryFunction(libraryfile,"vamp2_"+lc_functag);
      permfunc=s_loader->GetLibraryFunction(libraryfile,"set_permutation_"+lc_functag);
      if (ampfunc!=NULL && permfunc!=NULL) break;
    }
    if (ampfunc==NULL || permfunc==NULL) {
      PRINT_INFO("Didn't find functions");
      return NULL;
    }

    msg_Info()<<endl;
    PRINT_INFO("Initialising OpenLoopsReplace Virtual for "<<flavs<<": "<<lc_functag);
    return new OpenLoopsReplace_Virtual(pi, flavs, (Amp2Func) ampfunc,
                                 (PermutationFunc) permfunc, permutation, lc_functag);
  }
  else {
    if (OpenLoopsReplace_Interface::s_generate_list) {
      if (OpenLoopsReplace_Interface::s_shoppinglist.find(process)==OpenLoopsReplace_Interface::s_shoppinglist.end()) {
        ofstream list("OL_list.m", ios::app);
        list<<"(* "<<process<<" *)"<<endl;
        list<<"SubProcess["<<OpenLoopsReplace_Interface::s_shoppinglist.size()+1<<"] = {\n"
            <<" FeynArtsProcess -> {"
            <<OpenLoopsReplace_Interface::s_particles[map_flavs[0].HepEvt()].m_faname<<", "
            <<OpenLoopsReplace_Interface::s_particles[map_flavs[1].HepEvt()].m_faname<<"} -> {";
        for (size_t i=2; i<map_flavs.size()-1; ++i) {
          list<<OpenLoopsReplace_Interface::s_particles[map_flavs[i].HepEvt()].m_faname<<", ";
        }
        list<<OpenLoopsReplace_Interface::s_particles[map_flavs[map_flavs.size()-1].HepEvt()].m_faname<<"},\n"
            <<" SelectCoupling -> (Exponent[#, gQCD] == "<<pi.m_oqcd-1<<"+2*#2 &),\n"
            <<" SortExternal -> True,\n"
            <<" InsertFieldsOptions -> {Restrictions -> {ExcludeParticles -> {S[2|3], SV}, NoQuarkMixing}}";
        list<<"\n};\n"<<endl;
        list.close();
        PRINT_INFO("Generated list entry for "<<process);
        OpenLoopsReplace_Interface::s_shoppinglist.insert(process);
      }
      return new OpenLoopsReplace_Virtual(pi, flavs, dummyamp2func,
                                   dummypermfunc, permutation, "dummy");
    }
  }

  return NULL;
}