C++ Flavour::Mass方法代码示例

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.;
}

bool Soft_Cluster_Handler::
ClusterAnnihilation(Cluster * cluster,Flavour & had1,Flavour & had2) {
  int kfc1(int(cluster->GetTrip()->m_flav.Kfcode())); 
  int kfc2(int(cluster->GetAnti()->m_flav.Kfcode())); 
  kf_code kfc11(kfc1/1000),kfc12((kfc1-kfc11*1000)/100);
  kf_code kfc21(kfc2/1000),kfc22((kfc2-kfc21*1000)/100);
  Flavour fl1(kfc11), fl2(kfc12), fl3(kfc21), fl4(kfc22);
  fl1 = fl1.Bar();
  fl2 = fl2.Bar();
  Proto_Particle *pp1(new Proto_Particle(fl1,cluster->GetTrip()->m_mom/2.,'l'));
  Proto_Particle *pp2(new Proto_Particle(fl2,cluster->GetTrip()->m_mom/2.,'l'));
  Proto_Particle *pp3(new Proto_Particle(fl3,cluster->GetAnti()->m_mom/2.,'l'));
  Proto_Particle *pp4(new Proto_Particle(fl4,cluster->GetAnti()->m_mom/2.,'l'));
  bool order(ran->Get()>0.5?true:false);
  Cluster cluster1((order?pp3:pp4),pp1), cluster2((!order?pp4:pp3),pp2);
  Flavour_Pair pair1, pair2;
  pair1.first = cluster1.GetTrip()->m_flav;
  pair1.first = cluster1.GetAnti()->m_flav;
  pair2.first = cluster2.GetTrip()->m_flav;
  pair2.first = cluster2.GetAnti()->m_flav;
  double mass(cluster->Mass());
  double wt1(TransformWeight(&cluster1,had1,true,true));
  double wt2(TransformWeight(&cluster2,had2,true,true));
  if (wt1<0. || wt2<0.) return false;
  bool lighter1(wt1>0.), lighter2(wt2>0.);
  while (had1.Mass()+had2.Mass()>mass && lighter1 && lighter2) {
    lighter1 = wt1>0.;
    lighter2 = wt2>0.;
    if (wt1<=wt2) {
      if (lighter1)      wt1 = TransformWeight(&cluster1,had1,true,true);
      else if (lighter2) wt2 = TransformWeight(&cluster2,had2,true,true);
      else return false;
    }
    else {
      if (lighter2)      wt2 = TransformWeight(&cluster2,had2,true,true);
      else if (lighter1) wt1 = TransformWeight(&cluster1,had1,true,true);
      else return false;
    }
  }
  return true;
}

size_t Comix1to3::NHel(const Flavour& fl)
{
  switch(fl.IntSpin()) {
  case 0:
    return 1;
  case 1:
    return 2;
  case 2:
    if (IsZero(fl.Mass())) return 2;
    else return 3;
  default:
    THROW(not_implemented, "Comix not yet capable of spin > 1.");
    return 0;
  }
}

T3Channel::T3Channel(int _nin,int _nout,Flavour * fl,Flavour res) 
{  
  nin  = _nin; 
  nout = _nout;
  ms   = new double[nin+nout];
  for (int i=0;i<nin+nout;i++) ms[i] = ATOOLS::sqr(fl[i].Mass());
  rannum = 3*nout-4;
  rans   = new double[rannum];
  s      = smax  = pt2max = sqr(ATOOLS::rpa->gen.Ecms());
  pt2min = 0.0;
  E      = 0.5 * sqrt(s);
  name   = "T3-Channel";
  if (nout>2) name   = ToString(nin)+"->"+ToString(nout)+"_"+name;
  mass   = width = 0.; 
  type   = 0;
  if (res!=Flavour(kf_none)) {
    mass = res.Mass(); width = res.Width(); type = 1;
  }
  p_vegas = new Vegas(rannum,100,name,0);
}

Decay2Channel::Decay2Channel(int _nin,int _nout,const Flavour * fl,Flavour res) 
{  
  if (_nout != 2 || _nin!=1) {
    msg_Error()<<"Tried to initialize Decay2Channel with nout = "<<_nin<<" -> "<<_nout<<endl;
    abort();
  }
  nin  = _nin; nout = _nout;
  ms   = new double[nin+nout];
  for (short int i=0;i<nin+nout;i++) ms[i] = ATOOLS::sqr(fl[i].Mass());
  rannum = 2;
  rans   = new double[rannum];

  s      = smax  = pt2max = sqr(ATOOLS::rpa->gen.Ecms());
  pt2min = 0.;
  E      = 0.5 * sqrt(s);
  name   = "Decay2-Channel";
  mass   = width = 0.; 
  type   = 0;
  if (res!=Flavour(kf_none)) {
    mass = res.Mass(); width = res.Width(); type = 1;
  }
}

bool Hadron_Remnant::DecomposeHadron() 
{
  bool success=true;
  double Eb(p_beam->Energy());
  for (Particle_List::iterator pit=m_extracted.begin();
       pit!=m_extracted.end();++pit) {
    if ((*pit)->Momentum()[0]>Eb || (*pit)->Momentum()[0]<0.0) {
      msg_Error()<<"Hadron_Remnant::DecomposeHadron(): "
			 <<"Constituent energy out of range. \n   E_"
			 <<(*pit)->Flav()<<" = "<<(*pit)->Momentum()[0]
			 <<"."<<std::endl;
      success=false;
    }
    for (size_t j=0;j<m_constit.size();++j) {
      if ((*pit)->Flav()==m_constit[j]) {
	//std::cout<<METHOD<<" "<<success<<":"<<(*pit)->Flav()
	//	 <<" ("<<ValenceQuark(*pit)<<")"<<std::endl;
	if (success && ValenceQuark(*pit)) {
	  p_start = new Color_Dipole(*pit,&m_companions);  
	  p_start->Begin(ANTI((*pit)->Flav().IsAnti()))->
	    SetFlav(Opposite((*pit)->Flav()));
	  return success;
	}
      }
    }
  }
  Flavour    flav = m_constit[(size_t)(ran->Get()*3.)];
  Particle * part = new Particle(-1,flav); 
  part->SetStatus(part_status::active);
  part->SetFinalMass(flav.Mass());
  part->SetFlow(COLOR((qri::type)(flav.IsAnti())),Flow::Counter());
  //std::cout<<METHOD<<":"<<flav<<std::endl
  //	   <<"  "<<(*part)<<std::endl;
  p_start = new Color_Dipole(part,&m_companions);  
  p_start->Begin(ANTI(flav.IsAnti()))->SetFlav(Opposite(flav));
  m_companions.push_back(part);
  return success;
}

void Channel_Generator_Decays::GenerateMasses(int flag,Point** _plist,int pcount,
				       int& rannum,ofstream& sf)
{
  string * lm    = new string[pcount];
  string * momp  = new string[pcount];
  int    * sflag = new int[pcount];
  string mummy;
  string sum_s_i;
  string help;
  for (short int i=0;i<pcount;i++) {
    lm[i] = LinkedMasses(_plist[i]);
    mummy += lm[i];
    if (_plist[i]->left==0) {
      if (flag==0 || flag==10) AddToVariables(flag,lm[i],string("ms[")+GetMassIndex(lm[i])+string("]"),0,sf);
      //sf<<"  double s"<<lm[i]<<" = ms["<<lm[i]<<"];"<<endl;
      momp[i]  = string("p[") + GetMassIndex(lm[i]) + string("]");
      sflag[i] = 1;
      //sum_s_i  += string("-sqrt(s")+lm[i]+string(")");
      help    += lm[i];
    }
    else {
      CalcSmin(flag,"min",lm[i],sf,_plist[i]);
      momp[i]  = string("p") + Order(lm[i]);
      sflag[i] = 0;
    }
  }
  if (help.length()>0) {
    //CalcSmin(flag,"restmin",help,sf,0);
    //sum_s_i = string("-sqrt(s")+Order(help)+string("_restmin)");
    CalcSmin(flag,"min",help,sf,0);
    sum_s_i = string("-sqrt(s")+Order(help)+string("_min)");
  }
  int hit;
  double maxpole;
  double res;
  Flavour flav;
  string smax;
  for (;;) {
    hit = -1;
    maxpole = -1.;
    for (short int j=0;j<pcount;j++) {
      if (sflag[j]==0) {
	flav = _plist[j]->fl;
	res  = ATOOLS::sqr(flav.Width()*flav.Mass());
	if (!ATOOLS::IsZero(res) && Massive(flav)) {
	  if (1./res>maxpole) {
	    maxpole = 1./res;
	    hit = j;
	  }
	}
	else {
	  if (hit==-1) hit = j;
	}
      }
    }
    if (hit==-1) break;
    smax = string("sqr(sqrt(s")+Order(mummy)+string(")")+sum_s_i;
    
    for (short int j=0;j<pcount;j++) {
      if (sflag[j]==0 && j!=hit) {
	smax  += string("-sqrt(s")+Order(lm[j])+string("_min)");
      }
    }
    smax += string(")");
    
    AddToVariables(flag,lm[hit] +string("_max"),smax,0,sf);
    //sf<<"  double s"<<Order(lm[hit])<<"_max = "<<smax<<endl;

    int hi = 0;
    if (maxpole>0.) {
      hi = (_plist[hit]->fl).Kfcode();
      if (flag>=0) sf<<"  Flavour fl"<<lm[hit]<<" = "<<"Flavour((kf_code)("<<hi<<"));"<<endl;
    } 
    switch (flag) {
    case -11: case -1:
      if (maxpole>0.) {
	m_idc.push_back(string("MP")+ToString(hi)+string("_")+Order(lm[hit]));
      }
      else m_idc.push_back(string("MlP_")+Order(lm[hit]));
      break;
    case 0: case 10:
      sf<<"  Vec4D  "<<momp[hit]<<";"<<endl
	<<"  double s"<<Order(lm[hit])<<";"<<endl;
      if (maxpole>0.) {
	sf<<"  s"<<Order(lm[hit])
	  <<" = CE.MassivePropMomenta(fl"<<lm[hit]<<".Mass(),"<<"fl"<<lm[hit]<<".Width(),1,"
	  <<"s"<<Order(lm[hit])<<"_min,s"<<Order(lm[hit])<<"_max,ran["<<rannum<<"]);"<<endl;
      }
      else {
	sf<<"  s"<<Order(lm[hit])<<" = CE.MasslessPropMomenta(1.,s"<<Order(lm[hit])<<"_min,"
	  <<"s"<<Order(lm[hit])<<"_max,ran["<<rannum<<"]);"<<endl;
      }
      rannum++;
      break;
    default:
      string s(""); 
      for (int i=0;i<(int)lm[hit].length()-1;i++) s += string("p[")+GetMassIndex(lm[hit][i])+string("]+");
      s += string("p[")+GetMassIndex(lm[hit][lm[hit].length()-1])+string("]");
     
      AddToVariables(flag,lm[hit],s,1,sf);
//.........这里部分代码省略.........