C++ atools::Flavour类代码示例

bool Hard_Decay_Handler::Decays(const ATOOLS::Flavour& flav)
{
  if (flav.IsHadron()) return false;
  if (flav.Kfcode()==kf_tau && !m_decay_tau) return false;
  if (!flav.IsOn() || flav.IsStable()) return false;
  return true;
}

Primitive_Observable_Base *GetObservable(const Argument_Matrix &parameters)
{									
  if (parameters.size()<1) return NULL;
  if (parameters.size()==1) {
    if (parameters[0].size()<5) return NULL;
    int kf=ATOOLS::ToType<int>(parameters[0][0]);
    ATOOLS::Flavour flavour((kf_code)abs(kf));
    if (kf<0) flavour=flavour.Bar();
    std::string list=parameters[0].size()>5?parameters[0][5]:finalstate_list;
    return new Class(flavour,HistogramType(parameters[0][4]),
		     ATOOLS::ToType<double>(parameters[0][1]),
		     ATOOLS::ToType<double>(parameters[0][2]),
		     ATOOLS::ToType<int>(parameters[0][3]),list);
  }
  else if (parameters.size()<5) return NULL;
  double min=0.0, max=1.0;
  size_t bins=100;
  std::string list=finalstate_list, scale="Lin";
  ATOOLS::Flavour flavour;
  for (size_t i=0;i<parameters.size();++i) {
    if (parameters[i].size()<2) continue;
    if (parameters[i][0]=="FLAV") {
      int kf=ATOOLS::ToType<int>(parameters[i][1]);
      flavour=ATOOLS::Flavour((kf_code)abs(kf));
      if (kf<0) flavour=flavour.Bar();
    }
    else if (parameters[i][0]=="MIN") min=ATOOLS::ToType<double>(parameters[i][1]);
    else if (parameters[i][0]=="MAX") max=ATOOLS::ToType<double>(parameters[i][1]);
    else if (parameters[i][0]=="BINS") bins=ATOOLS::ToType<int>(parameters[i][1]);
    else if (parameters[i][0]=="SCALE") scale=parameters[i][1];
    else if (parameters[i][0]=="LIST") list=parameters[i][1];
  }
  return new Class(flavour,HistogramType(scale),min,max,bins,list);
}									

Pol_Info::Pol_Info(const ATOOLS::Flavour& fl)
{
  int dof = 1;
  pol_type='s';
  if(fl.IsFermion())                                 { dof = 2;pol_type='h';};
  if(fl.IsVector() &&  ATOOLS::IsZero(fl.Mass()))  { dof = 2;pol_type='c';}
  if(fl.IsVector() && !ATOOLS::IsZero(fl.Mass()))  {

#ifdef Explicit_Pols
    dof=3;
#else
    dof=1;
#endif
    pol_type='c';
  }
  if(fl.IsTensor()) dof=5;
  Init(dof);

  if(!fl.IsTensor()) {
    int tf[3]  = {mt::p_m, mt::p_p, mt::p_l };
    for(int j=0;j<dof;j++){
      type[j]   = tf[j];
      factor[j] = 1.;
    }
  } 
  else {
    type[0]=mt::p_t1;factor[0]=1.;
    type[1]=mt::p_t2;factor[1]=1.;
    type[2]=mt::p_t3;factor[2]=1.;
    type[3]=mt::p_t4;factor[3]=1.;
    type[4]=mt::p_t5;factor[4]=1.;
  }
}

double Kinematics_FF::GetKT2(const double &Q2,const double &y,const double &z,
			     const double &mi2,const double &mj2,const double &mk2,
			     const ATOOLS::Flavour &fla,const ATOOLS::Flavour &flc) const
{
  double pipj=(Q2-mi2-mj2-mk2)*y;
  if (m_evolscheme==0) {
    return pipj*z*(1.0-z)-sqr(1.0-z)*mi2-sqr(z)*mj2;
  }
  double kt2=pipj*z*(1.0-z);
  if (fla.IsFermion()) kt2=pipj*(flc.IsVector()?(1.0-z):z);
  else if (flc.IsFermion()) kt2=pipj;
  return kt2;
}

bool CF_EW_FFZ::AllowSpec(const ATOOLS::Flavour &fl)
{
    if (m_cfl.IntCharge()==0) return fl.Charge();

    switch (m_type) {
    case cstp::FF:
    case cstp::II:
        return fl.IntCharge()*m_cfl.IntCharge()<0;
    case cstp::FI:
    case cstp::IF:
        return fl.IntCharge()*m_cfl.IntCharge()>0;
    default:
        return false;
    }
}

bool CF_QED::AllowSpec(const ATOOLS::Flavour &fl) 
{
  if (!fl.Strong() && fl.Mass()>10.0) return false;
  if (m_cfl.IntCharge()==0) return fl.Charge();
  
  switch (m_type) {
  case cstp::FF:
  case cstp::II:
    return fl.IntCharge()*m_cfl.IntCharge()<0;
  case cstp::FI:
  case cstp::IF:
    return fl.IntCharge()*m_cfl.IntCharge()>0;
  default:
    return false;
  }
}

double Kinematics_FF::GetY(const double &Q2,const double &kt2,const double &z,
			   const double &mi2,const double &mj2,const double &mk2,
			   const ATOOLS::Flavour &fla,const ATOOLS::Flavour &flc,
			   const bool force) const
{
  if (!force && (z<=0.0 || z>=1.0 || Q2<=mi2+mj2+mk2)) return -1.0;
  if (m_evolscheme==0) {
    return (kt2/(z*(1.0-z))+(1.0-z)/z*mi2+z/(1.0-z)*mj2)/(Q2-mi2-mj2-mk2);
  }
  if (fla.IsFermion()) {
    if (flc.IsFermion()) return kt2/z/(Q2-mi2-mj2-mk2);
    return kt2/(1.0-z)/(Q2-mi2-mj2-mk2);
  }
  if (flc.IsFermion()) return kt2/(Q2-mi2-mj2-mk2);
  return kt2/(z*(1.0-z))/(Q2-mi2-mj2-mk2);
}

double Kinematics_II::GetKT2(const double &Q2,const double &y,const double &z,
			     const double &ma2,const double &mi2,const double &mb2,
			     const ATOOLS::Flavour &flb,const ATOOLS::Flavour &flc) const
{
  double pipj=(Q2-ma2-mi2-mb2)*y/z;
  if (m_evolscheme==0) {
    return pipj*(1.0-z)-mi2-sqr(1.0-z)*ma2;
  }
  double kt2=pipj*(1.0-z);
  if (flc.IsFermion()) kt2=pipj;
  return kt2;
}

double Hadron_Remnant::GetXPDF(ATOOLS::Flavour flavour,double scale) 
{
  double cut, x;
  cut=2.0*(flavour.HadMass()+m_hardpt.PPerp()/
	   sqr(m_companions.size()))/p_beam->OutMomentum()[0];
  // assume heavy flavours have been pair-produced
  // => scale should be approximately (2m)^2
  scale=Max(scale,4.0*sqr(flavour.Mass()));
  if (scale<p_pdfbase->Q2Min()) {
    msg_Error()<<"Hadron_Remnant::GetXPDF("<<flavour<<","<<scale<<"): "
		       <<"Scale under-runs minimum given by PDF: "
		       <<scale<<" < "<<p_pdfbase->Q2Min()<<std::endl;
    scale=1.001*p_pdfbase->Q2Min();
  } 
  unsigned int xtrials, pdftrials=0;
  while (true) {
    ++pdftrials;
    xtrials=0;
    do { 
      ++xtrials;
      x=m_xrem*ran->Get();
      if (xtrials>=m_maxtrials) {
	x=Min(cut,0.999999*p_pdfbase->RescaleFactor());
	break;
      }
      if (x>p_pdfbase->RescaleFactor()) continue;
    } while (x<cut);	
    if (x>p_pdfbase->XMin() && x<p_pdfbase->XMax()) {
      p_pdfbase->Calculate(x,scale);
    }
    else {
      m_xscheme=0; return 0.01;
    }
    if (pdftrials>=m_maxtrials) { m_xscheme=0; return 0.01; }
    if (p_pdfbase->GetXPDF(flavour)/x>ran->Get()) return x;
  } 
  return 0.0;
}

Channel_Interface::Channel_Interface(int nin,int nout,ATOOLS::Flavour *flavour,ATOOLS::Flavour res):
  Single_Channel(nin,nout,flavour)
{  
  if (nout != 2 || nin!=2) {
    msg_Error()<<"Channel_Interface::Channel_Interface(..): "
		       <<"Cannot handle "<<nin<<" -> "<<nout<<" processes. Abort."<<std::endl;
    exit(169);
  }
  nin  = nin; nout = nout;
  ms = new double[nin+nout];
  for (short int i=0;i<nin+nout;i++) ms[i] = ATOOLS::sqr(flavour[i].Mass());
  rannum = 3;
  rans = new double[rannum];
  s = smax = pt2max = ATOOLS::sqr(ATOOLS::rpa->gen.Ecms());
  pt2min = 0.;
  E = 0.5*sqrt(s);
  name = "Channel Interface";
  mass = width = 0.; 
  type = 0;
  if (res!=ATOOLS::Flavour(kf_none)) {
    mass = res.Mass(); width = res.Width(); type = 1;
  }
}

ATOOLS::Flavour Hadron_Remnant::Opposite(ATOOLS::Flavour flav) const
{
  bool found=false;
  kf_code rem[2];
  for (short unsigned int i=0,j=0;i<3;++i) {
    if (m_constit[i]==flav && !found) found=true;
    else rem[j++]=m_constit[i].Kfcode();
  }
  Flavour anti=Flavour((kf_code)(abs(rem[0])*1000+abs(rem[1])*100+3));
  if (rem[0]!=rem[1]) {
    if (ran->Get()<0.25) 
      anti=Flavour((kf_code)(abs(rem[0])*1000+abs(rem[1])*100+1));
  }
  else {
    anti=Flavour((kf_code)(abs(rem[0])*1100+3));
  }
  if (flav.IsAnti()) anti=anti.Bar();
  return anti;
}

 double GetXPDF(const ATOOLS::Flavour& infl) {
   if (m_x>m_xmax || m_rescale<0.) return 0.;
   if (!(m_x>=0.0 && m_x<=1.0)) {
     PRINT_INFO("PDF called with x="<<m_x);
     return 0.;
   }
   int cteqindex;
   switch (infl.Kfcode()) {
   case kf_gluon: cteqindex=0;                  break;
   case kf_d:     cteqindex=m_anti*int(infl)*2; break;
   case kf_u:     cteqindex=m_anti*int(infl)/2; break;
   default:                cteqindex=m_anti*int(infl);   break;
   }
   if (!m_calculated[5-cteqindex]) {
     m_f[5-cteqindex]=ct10pdf_(cteqindex,m_x,m_Q)*m_x; 
     m_calculated[5-cteqindex]=true;
   }
   return m_rescale*m_f[5-cteqindex];     
 }

double Hadron_Remnant::MinimalEnergy(const ATOOLS::Flavour &flavour) 
{
  if (!m_initialized) {
    if (!flavour.Strong()) {
      return p_beam->Beam().HadMass();
    }
    bool found(false);
    kf_code di[3];
    if (flavour.IsQuark()) {
      short unsigned int j=0;
      for (Flavour_Vector::const_iterator flit(m_constit.begin());
	   flit!=m_constit.end();++flit) {
	if (found || flavour!=*flit) di[j++]=flit->Kfcode();
	else found=true;
      }
    }
    Flavour difl;
    if (!found || flavour.IsGluon()) {
      int single=-1;
      for (size_t i=0;i<m_constit.size();++i) {
	for (size_t j=i+1;j<m_constit.size();++j) {
	  if (m_constit[i]==m_constit[j]) { 
	    single=j; 
	    break; 
	  }
	}
	if (single>0) break;
      }
      Flavour fl(m_constit[single]);
      for (short unsigned int j=0, i=0;i<3;i++) 
	if (i!=single) di[j++]=m_constit[i].Kfcode();
      if (di[0]>di[1]) difl=Flavour((kf_code)(di[0]*1000+di[1]*100+1));
      else if (di[1]>di[0]) difl=Flavour((kf_code)(di[1]*1000+di[0]*100+1));
      else difl=Flavour((kf_code)(di[0]*1100+3));
      if (m_constit[single].IsAnti()) difl=difl.Bar();
      return difl.Mass()+fl.Mass()+flavour.Bar().Mass();
    }
    if (di[0]>di[1]) difl=Flavour((kf_code)(di[0]*1000+di[1]*100+1));
    else if (di[1]>di[0]) difl=Flavour((kf_code)(di[1]*1000+di[0]*100+1));
    else difl=Flavour((kf_code)(di[0]*1100+3));
    if (m_constit[0].IsAnti()) difl=difl.Bar();
    return difl.Mass();
  }
  else {
    if (flavour.IsQuark()) return flavour.Bar().Mass();
  }
  return 0.;
}

double MI_Handler::Mass(const ATOOLS::Flavour &fl) const
{
  return fl.Mass();
}

double Hard_Decay_Handler::Mass(const ATOOLS::Flavour &fl) const
{
  if (m_usemass==0) return fl.Mass();
  if (m_decmass.find(fl)!=m_decmass.end()) return fl.Mass(true);
  return fl.Mass();
}