Ticket #223: C7_615.C

#include "PHASIC++/Channels/Single_Channel.H"
#include "ATOOLS/Org/Run_Parameter.H"
#include "PHASIC++/Channels/Channel_Elements.H"
#include "PHASIC++/Channels/Vegas.H"

using namespace PHASIC;
using namespace ATOOLS;

namespace PHASIC {
  class C7_615 : public Single_Channel {
    double m_amct,m_alpha,m_ctmax,m_ctmin;
    Info_Key m_kTC23_0_378_1_26_5_4,m_kTC24_0_8_1_6_7_2345,m_kTC24_0_8_1__23457_6,m_kTC_0_78_1_26_3_45,m_kTC_0_78_1_6_345_2,m_kTC_0__1__8_234567,m_kZS_0;
    Vegas* p_vegas;
  public:
    C7_615(int,int,Flavour*,Integration_Info * const);
    ~C7_615();
    void   GenerateWeight(Vec4D *,Cut_Data *);
    void   GeneratePoint(Vec4D *,Cut_Data *,double *);
    void   AddPoint(double);
    void   MPISync()                 { p_vegas->MPISync(); }
    void   Optimize()                { p_vegas->Optimize(); } 
    void   EndOptimize()             { p_vegas->EndOptimize(); } 
    void   WriteOut(std::string pId) { p_vegas->WriteOut(pId); } 
    void   ReadIn(std::string pId)   { p_vegas->ReadIn(pId); } 
    void   ISRInfo(int &,double &,double &);
    std::string ChID();
  };
}

extern "C" Single_Channel * Getter_C7_615(int nin,int nout,Flavour* fl,Integration_Info * const info) {
  return new C7_615(nin,nout,fl,info);
}

void C7_615::GeneratePoint(Vec4D * p,Cut_Data * cuts,double * _ran)
{
  double *ran = p_vegas->GeneratePoint(_ran);
  for(int i=0;i<rannum;i++) rans[i]=ran[i];
  Vec4D p2345678=p[0]+p[1];
  double s2345678_max = p2345678.Abs2();
  double s234567_max = sqr(sqrt(s2345678_max)-sqrt(ms[8]));
  double s6 = ms[6];
  double s68_min = cuts->Getscut(std::string("68"));
  double s23457_max = sqr(sqrt(s2345678_max)-sqrt(s68_min));
  double s678_min = cuts->Getscut(std::string("678"));
  double s2345_max = sqr(sqrt(s2345678_max)-sqrt(s678_min));
  double s2 = ms[2];
  double s2678_min = cuts->Getscut(std::string("2678"));
  double s345_max = sqr(sqrt(s2345678_max)-sqrt(s2678_min));
  double s23678_min = cuts->Getscut(std::string("23678"));
  double s45_max = sqr(sqrt(s2345678_max)-sqrt(s23678_min));
  double s4 = ms[4];
  double s5 = ms[5];
  double s45_min = cuts->Getscut(std::string("45"));
  Vec4D  p45;
  double s45 = CE.ThresholdMomenta(1.5,4.*sqrt(s45_min),s45_min,s45_max,ran[0]);
  double s3 = ms[3];
  double s345_min = cuts->Getscut(std::string("345"));
  s345_min = Max(s345_min,sqr(sqrt(s3)+sqrt(s45)));
  Vec4D  p345;
  double s345 = CE.ThresholdMomenta(1.5,2.*sqrt(s345_min),s345_min,s345_max,ran[1]);
  double s2345_min = cuts->Getscut(std::string("2345"));
  s2345_min = Max(s2345_min,sqr(sqrt(s345)+sqrt(s2)));
  Vec4D  p2345;
  double s2345 = CE.ThresholdMomenta(1.5,2.*sqrt(s2345_min),s2345_min,s2345_max,ran[2]);
  double s7 = ms[7];
  double s23457_min = cuts->Getscut(std::string("23457"));
  s23457_min = Max(s23457_min,sqr(sqrt(s7)+sqrt(s2345)));
  Vec4D  p23457;
  double s23457 = CE.ThresholdMomenta(1.5,2.*sqrt(s23457_min),s23457_min,s23457_max,ran[3]);
  double s234567_min = cuts->Getscut(std::string("234567"));
  s234567_min = Max(s234567_min,sqr(sqrt(s23457)+sqrt(s6)));
  Vec4D  p234567;
  double s234567 = CE.ThresholdMomenta(1.5,2.*sqrt(s234567_min),s234567_min,s234567_max,ran[4]);
  double s8 = ms[8];
  m_ctmax = cuts->cosmax[0][8];
  m_ctmin = cuts->cosmin[0][8];
  CE.TChannelMomenta(p[0],p[1],p[8],p234567,s8,s234567,0.,m_alpha,m_ctmax,m_ctmin,m_amct,0,ran[5],ran[6]);
  double tmass206 = Flavour((kf_code)(24)).Mass();
  Vec4D  p0_8 = p[0]-p[8];
  CE.TChannelMomenta(p0_8,p[1],p23457,p[6],s23457,s6,tmass206,m_alpha,1.,-1.,m_amct,0,ran[7],ran[8]);
  double tmass202 = Flavour((kf_code)(24)).Mass();
  Vec4D  p1_6 = p[1]-p[6];
  CE.TChannelMomenta(p0_8,p1_6,p[7],p2345,s7,s2345,tmass202,m_alpha,1.,-1.,m_amct,0,ran[9],ran[10]);
  Vec4D  p0_78 = p0_8-p[7];
  CE.TChannelMomenta(p0_78,p1_6,p345,p[2],s345,s2,0.,m_alpha,1.,-1.,m_amct,0,ran[11],ran[12]);
  Vec4D  p1_26 = p1_6-p[2];
  CE.TChannelMomenta(p0_78,p1_26,p[3],p45,s3,s45,0.,m_alpha,1.,-1.,m_amct,0,ran[13],ran[14]);
  double tmass204 = Flavour((kf_code)(23)).Mass();
  Vec4D  p0_378 = p0_87-p[3];
  CE.TChannelMomenta(p0_378,p1_26,p[5],p[4],s5,s4,tmass204,m_alpha,1.,-1.,m_amct,0,ran[15],ran[16]);
}

void C7_615::GenerateWeight(Vec4D* p,Cut_Data * cuts)
{
  double wt = 1.;
  Vec4D p2345678=p[0]+p[1];
  double s2345678_max = p2345678.Abs2();
  double s234567_max = sqr(sqrt(s2345678_max)-sqrt(ms[8]));
  double s6 = ms[6];
  double s68_min = cuts->Getscut(std::string("68"));
  double s23457_max = sqr(sqrt(s2345678_max)-sqrt(s68_min));
  double s678_min = cuts->Getscut(std::string("678"));
  double s2345_max = sqr(sqrt(s2345678_max)-sqrt(s678_min));
  double s2 = ms[2];
  double s2678_min = cuts->Getscut(std::string("2678"));
  double s345_max = sqr(sqrt(s2345678_max)-sqrt(s2678_min));
  double s23678_min = cuts->Getscut(std::string("23678"));
  double s45_max = sqr(sqrt(s2345678_max)-sqrt(s23678_min));
  double s4 = ms[4];
  double s5 = ms[5];
  double s45_min = cuts->Getscut(std::string("45"));
  Vec4D  p45 = p[4]+p[5];
  double s45 = dabs(p45.Abs2());
  wt *= CE.ThresholdWeight(1.5,4.*sqrt(s45_min),s45_min,s45_max,s45,rans[0]);
  double s3 = ms[3];
  double s345_min = cuts->Getscut(std::string("345"));
  s345_min = Max(s345_min,sqr(sqrt(s3)+sqrt(s45)));
  Vec4D  p345 = p[3]+p[4]+p[5];
  double s345 = dabs(p345.Abs2());
  wt *= CE.ThresholdWeight(1.5,2.*sqrt(s345_min),s345_min,s345_max,s345,rans[1]);
  double s2345_min = cuts->Getscut(std::string("2345"));
  s2345_min = Max(s2345_min,sqr(sqrt(s345)+sqrt(s2)));
  Vec4D  p2345 = p[2]+p[3]+p[4]+p[5];
  double s2345 = dabs(p2345.Abs2());
  wt *= CE.ThresholdWeight(1.5,2.*sqrt(s2345_min),s2345_min,s2345_max,s2345,rans[2]);
  double s7 = ms[7];
  double s23457_min = cuts->Getscut(std::string("23457"));
  s23457_min = Max(s23457_min,sqr(sqrt(s7)+sqrt(s2345)));
  Vec4D  p23457 = p[2]+p[3]+p[4]+p[5]+p[7];
  double s23457 = dabs(p23457.Abs2());
  wt *= CE.ThresholdWeight(1.5,2.*sqrt(s23457_min),s23457_min,s23457_max,s23457,rans[3]);
  double s234567_min = cuts->Getscut(std::string("234567"));
  s234567_min = Max(s234567_min,sqr(sqrt(s23457)+sqrt(s6)));
  Vec4D  p234567 = p[2]+p[3]+p[4]+p[5]+p[6]+p[7];
  double s234567 = dabs(p234567.Abs2());
  wt *= CE.ThresholdWeight(1.5,2.*sqrt(s234567_min),s234567_min,s234567_max,s234567,rans[4]);
  double s8 = ms[8];
  m_ctmax = cuts->cosmax[0][8];
  m_ctmin = cuts->cosmin[0][8];
  if (m_kTC_0__1__8_234567.Weight()==ATOOLS::UNDEFINED_WEIGHT)
    m_kTC_0__1__8_234567<<CE.TChannelWeight(p[0],p[1],p[8],p234567,0.,m_alpha,m_ctmax,m_ctmin,m_amct,0,m_kTC_0__1__8_234567[0],m_kTC_0__1__8_234567[1]);
  wt *= m_kTC_0__1__8_234567.Weight();

  rans[5]= m_kTC_0__1__8_234567[0];
  rans[6]= m_kTC_0__1__8_234567[1];
  double tmass206 = Flavour((kf_code)(24)).Mass();
  Vec4D  p0_8 = p[0]-p[8];
  if (m_kTC24_0_8_1__23457_6.Weight()==ATOOLS::UNDEFINED_WEIGHT)
    m_kTC24_0_8_1__23457_6<<CE.TChannelWeight(p0_8,p[1],p23457,p[6],tmass206,m_alpha,1.,-1.,m_amct,0,m_kTC24_0_8_1__23457_6[0],m_kTC24_0_8_1__23457_6[1]);
  wt *= m_kTC24_0_8_1__23457_6.Weight();

  rans[7]= m_kTC24_0_8_1__23457_6[0];
  rans[8]= m_kTC24_0_8_1__23457_6[1];
  double tmass202 = Flavour((kf_code)(24)).Mass();
  Vec4D  p1_6 = p[1]-p[6];
  if (m_kTC24_0_8_1_6_7_2345.Weight()==ATOOLS::UNDEFINED_WEIGHT)
    m_kTC24_0_8_1_6_7_2345<<CE.TChannelWeight(p0_8,p1_6,p[7],p2345,tmass202,m_alpha,1.,-1.,m_amct,0,m_kTC24_0_8_1_6_7_2345[0],m_kTC24_0_8_1_6_7_2345[1]);
  wt *= m_kTC24_0_8_1_6_7_2345.Weight();

  rans[9]= m_kTC24_0_8_1_6_7_2345[0];
  rans[10]= m_kTC24_0_8_1_6_7_2345[1];
  Vec4D  p0_78 = p0_8-p[7];
  if (m_kTC_0_78_1_6_345_2.Weight()==ATOOLS::UNDEFINED_WEIGHT)
    m_kTC_0_78_1_6_345_2<<CE.TChannelWeight(p0_78,p1_6,p345,p[2],0.,m_alpha,1.,-1.,m_amct,0,m_kTC_0_78_1_6_345_2[0],m_kTC_0_78_1_6_345_2[1]);
  wt *= m_kTC_0_78_1_6_345_2.Weight();

  rans[11]= m_kTC_0_78_1_6_345_2[0];
  rans[12]= m_kTC_0_78_1_6_345_2[1];
  Vec4D  p1_26 = p1_6-p[2];
  if (m_kTC_0_78_1_26_3_45.Weight()==ATOOLS::UNDEFINED_WEIGHT)
    m_kTC_0_78_1_26_3_45<<CE.TChannelWeight(p0_78,p1_26,p[3],p45,0.,m_alpha,1.,-1.,m_amct,0,m_kTC_0_78_1_26_3_45[0],m_kTC_0_78_1_26_3_45[1]);
  wt *= m_kTC_0_78_1_26_3_45.Weight();

  rans[13]= m_kTC_0_78_1_26_3_45[0];
  rans[14]= m_kTC_0_78_1_26_3_45[1];
  double tmass204 = Flavour((kf_code)(23)).Mass();
  Vec4D  p0_378 = p0_87-p[3];
  if (m_kTC23_0_378_1_26_5_4.Weight()==ATOOLS::UNDEFINED_WEIGHT)
    m_kTC23_0_378_1_26_5_4<<CE.TChannelWeight(p0_378,p1_26,p[5],p[4],tmass204,m_alpha,1.,-1.,m_amct,0,m_kTC23_0_378_1_26_5_4[0],m_kTC23_0_378_1_26_5_4[1]);
  wt *= m_kTC23_0_378_1_26_5_4.Weight();

  rans[15]= m_kTC23_0_378_1_26_5_4[0];
  rans[16]= m_kTC23_0_378_1_26_5_4[1];
  double vw = p_vegas->GenerateWeight(rans);
  if (wt!=0.) wt = vw/wt/pow(2.*M_PI,7*3.-4.);

  weight = wt;
}

C7_615::C7_615(int nin,int nout,Flavour* fl,Integration_Info * const info)
       : Single_Channel(nin,nout,fl)
{
  name = std::string("C7_615");
  rannum = 17;
  rans  = new double[rannum];
  m_amct  = 1.;
  m_alpha = .9;
  m_ctmax = 1.;
  m_ctmin = -1.;
  m_kTC23_0_378_1_26_5_4.Assign(std::string("TC23_0_378_1_26_5_4"),2,0,info);
  m_kTC24_0_8_1_6_7_2345.Assign(std::string("TC24_0_8_1_6_7_2345"),2,0,info);
  m_kTC24_0_8_1__23457_6.Assign(std::string("TC24_0_8_1__23457_6"),2,0,info);
  m_kTC_0_78_1_26_3_45.Assign(std::string("TC_0_78_1_26_3_45"),2,0,info);
  m_kTC_0_78_1_6_345_2.Assign(std::string("TC_0_78_1_6_345_2"),2,0,info);
  m_kTC_0__1__8_234567.Assign(std::string("TC_0__1__8_234567"),2,0,info);
  m_kZS_0.Assign(std::string("ZS_0"),2,0,info);
  p_vegas = new Vegas(rannum,100,name);
}

C7_615::~C7_615()
{
  delete p_vegas;
}

void C7_615::ISRInfo(int & type,double & mass,double & width)
{
  type  = 2;
  mass  = 0;
  width = 0.;
}

void C7_615::AddPoint(double Value)
{
  Single_Channel::AddPoint(Value);
  p_vegas->AddPoint(Value,rans);
}
std::string C7_615::ChID()
{
  return std::string("CG2$MTH_2345$MTH_234567$MTH_23457$MTH_345$MTH_45$TC23_0_378_1_26_5_4$TC24_0_8_1_6_7_2345$TC24_0_8_1__23457_6$TC_0_78_1_26_3_45$TC_0_78_1_6_345_2$TC_0__1__8_234567$ZS_0$");
}